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Magnetogasdvnamics of Hypersonic 
Couette Flow 


Z. O. BLEVISS* 
Douglas Aircraft Company, Inc. 


SUMMARY 


When the gas in a high-speed laminar boundary layer becomes 
jonized due to viscous heating, it becomes an electrical conductor 
aud it can be made to interact with a magnetic field to alter skin 
friction, heat transfer, total drag, ete. The calculation of mag- 
netogasdynamic effects on high-speed laminar boundary layers 
is extremely complex and no solutions with realistic assumptions 
have been obtained. A very simple pure shear flow, usually 
called Couette flow, is investigated in this paper because it con- 
tains many important features of boundary-layer flow and_ be- 
cause it can be solved exactly with a minimum of assumptions 
about the gas. 

Couette flow is the fluid motion produced when a viscous fluid 
contained between two infinitely extended plane parallel walls is 
set into steady motion by the relative movement of the walls in 
their own planes. The exact solution for compressible Couette 
flow but with no magnetic field was given by Illingworth.? 
Liepmann and Bleviss* showed how easily effects like equilibrium 
dissociation and ionization can be included in the solution and 
drew attention to the fact that a study of this type of flow is very 
helpful in understanding the flow in a high-speed laminar bound- 
ary layer. 

In the present paper, the Couette flow problem is extended to 
include a uniform externally imposed magnetic field normal to 
the walls. Assuming thermodynamic equilibrium and reason- 
able variations of electrical conductivity, viscosity, and Prandtl 
Number with temperature, numerical solutions are presented for 
air for the case of the insulated wall (no heat transfer to the wall) 
for Mach Numbers from 10 to 30 and for the case of heat trans- 
fer for Mach Numbers of 20 and 30. The effects of the magnetic 
field upon the velocity, temperature, electrical current density, 
and induced magnetic field distributions and upon the skin fric- 
tion, heat transfer, and total drag—i.e., skin friction plus the 
Magnetic stress due to the interaction of the electrically conduct- 
ing fluid with the magnetic field—are shown. 


Received March 7, 1958. Revised and received June 25, 1958. 
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Comparing results with magnetic field to results without mag- 
netic field for a given high-speed laminar boundary layer, the 
present investigation shows that relatively weak magnetic fields 
produce large increases in the total drag, large reductions in the 
skin friction, and at the same time have relatively little effect on 
the heat transfer. Whereas the total drag without magnetic 
field is skin-friction drag, the total drag with magnetic field is 
primarily magnetic drag. This implies that if the comparison 
is made on the basis of the same total drag the heat transfer for 
the magnetic case will be much less than that for the nonmagnetic 
case. 

The expected qualitative and quantitative applicability of the 
Couette flow results to the boundary layer is discussed. 


SYMBOLS 


magnetic field strength (magnetic flux density) 
m.k.s. unit is weber/meter?. Unit used in 
presentation of results is gauss (1 weber + 
meter? = 104 gauss). 

imposed magnetic field (see Fig. 1) 

Cartesian components of B (see Fig. 1) 

t./{(1/2)p,, U?] = total drag coefficient 

7,/{1/2)p,, U2] = skin-friction coefficient 

gy /(h, h,) p..U = heat-transfer coefficient 

Ga/(h, — h,) p.U = heat-transfer coefficient 
for moving wall 

specific heat at constant pressure 

electric field intensity (volt /meter) 

specific enthalpy (enthalpy per unit mass 

h/he« 


i. maximum enthalpy for zero magnetic field 


J current density (ampere/meter?) 
J = total current per unit length in x direction 
(ampere /meter) 
thermal conductivity 
free-stream Mach Number 
V Bod = Hartmann Number 
fluid mechanical pressure 
fluid mechanical stress tensor 
components of fluid mechanical stress tensor 
[see Eq. (Gc)] 
nC,/k = Prandtl Number 
heat flux 
woo 6U = magnetic Reynolds Number 


max 
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Fic. 1. Couette flow with transverse magnetic field. 

R; = poli/no = fluid Reynolds Number based on 
free-stream properties and 6 

= absolute temperature 

D/T. 

Z = Maxwell stress tensor 

Tix. = components of Maxwell stress tensor [see Eq. 
(6d)] 

u = velocity component in x direction (see Fig. 1) 

= u/U 

U = velocity of moving wall or free stream (see 
Fig. 1) 

V = total velocity 

x, Y, 2 = right-hand Cartesian coordinate system (see 
Fig. 1) 

= y/é 

Y = ratio of specific heats in free stream = 1.40 

6 = distance between walls or boundary-layer thick- 
ness (see Fig. 1) 

bi, = components of unit tensor (Kronecker delta) 

€ = dielectric constant (farad/meter) 

n = viscosity 

n/N 

= magnetic permeability (hency/meter) 

p = density 

a = electrical conductivity (mho/meter) 

em = maximum electrical conductivity for zero mag- 
netic field 

a = «/oomar 

T = shear stress in Couette flow 

7 = 

“ = viscous part of fluid mechanical stress tensor 

Tik = components of viscous part of fluid mechanical 
stress tensor 

Subscripts 

( )o = zero magnetic field 

= recovery conditions—i.e., conditions at station- 
ary wall with zero heat transfer 

(> = conditions at stationary wall or body (see Fig. 1) 

om = conditions at moving wall or in free stream 
a (see Fig. 1) 
= avector 
INTRODUCTION 


N HYPERSONIC SPEEDS, the gas in a high-speed 
laminar boundary layer becomes ionized due to 
viscous heating. The gas then becomes an electrical 
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conductor and it can be made to interact with a mag 


netic field to alter the skin friction, heat transfer, toty} 


drag, etc. The exact calculation of magnetogasdynam; 
effects on high-speed laminar boundary layers is ¢. 
tremely complicated and no solutions with realistic 
assumptions have been obtained. Some indication 9j 
these difficulties can be found in reference 1, where th 
incompressible, laminar, flat-plate boundary | 
with a uniform magnetic field normal to the plate is 
treated approximately. Because of these difficulties 
a very simple pure shear flow that contains many oj 
the important features of boundary-layer flow is iy. 
vestigated here. The really significant advantag 
gained by resorting to this simpler problem is that jt 


can be solved exactly for all Mach and Reynolds Num. 
bers with a minimum number of assumptions about the 
gas. This simple flow, usually called Couette flow, then | 
represents an excellent test case to study the effects oj. 
a magnetic field on recovery temperature, skin friction, 
heat transfer, etc. because the physics of the problem | 
is not hidden behind a formidable mathematical and 
numerical apparatus. 

Couette flow is the fluid motion produced when a | 
viscous fluid contained between two infinitely extended 
plane parallel walls is set into steady motion by the 
relative movement of the walls in their own planes (see 
lower sketch in Fig. 1). Consider the two walls at 
y = 0 and y = 6 and let the upper wall move with 
constant velocity LU’ relative to the lower wall which | 
will be taken to be stationary. In the comparison | 
with the boundary layer the moving wall corresponds | 
to the free stream outside the boundary layer, the sta- 
tionary wall corresponds to the body, and 6 is the 
boundary-layer thickness. When the electromagnetic 
boundary conditions are considered it will be found | 
convenient to interpret Couette flow as an approxima. | 
tion to the fluid motion produced in a viscous fluid con- | 
tained in the annular region between two coaxial circu- 
lar* cylinders of infinite length sliding relative to each | 
other along their common axis (see upper sketch in 
Fig. 1), the thickness 6 of the annular region being small | 
compared with the cylinder radii. This configuration | 
corresponds to a boundary layer on the slender after: | 
body of a body of revolution. The notation in this 
report has been chosen to correspond as closely as 
possible to the usual boundary-layer notation. The 
stationary wall will be referred to as the ‘‘wall’’ and all 
quantities there will be denoted by the subscript w—i-. 
Tw, Ty, etc. Except where specifically stated other 
wise, the upper wall (y 5) will correspond to the 
“free stream’ and will be denoted by the subscript 
© (the ‘‘free-stream’”’ velocity is denoted by U). 

The exact solution for compressible Couette flow 
when no electromagnetic effects are considered was 
given by Illingworth.2. Liepmann and Bleviss* showed 


* The cylinder cross sections need not be circular. Howevel, 
the normal distance 6 between the cylinders must be constail 
around the circumference and must be small compared with the 
smallest radius of curvature. For convenience, the cylinder 
cross sections will always be taken circular. 


how 
joniza' 
drew 
flow 
high-s 

In t 
extend 
netic fi 
eylind 
Unless 
tached 
most 11 


Bef« 
Couett 
of the 
feature 
flow of 
electric 
The ra 
is used 

Fror 


Ohm 


where 
The 
the mo 
netic fie 
The | 
The co: 


A conv 


The dy 
the secc 


P is th 
compon 


where ; 
referenc 
tensor \ 


¥ 602 
| 
\ 
\ i} 
>U 
uly) 
8 y 
x 


Mag 
total 
Namic 
is 
alistic 
ion 
re the 


layer 


ate js 
ulties 
ny oj 
iS in- 
Ntage 
hat it 
Num- 


it the | 


, then 
cts of 
ction, 
and 


len 
nded 
y the 
(see 
Is at 
with 


vhich 
rison 


onds 
sta- 
> the 
netic 


ound | 


ima- | 


con- 


ircu- 


each 
h in 


mal! | 


ition 


fter- | 


this 
y as 
The 
d all 
4.€., 


her- 


MAGNETOGASDYNAMICS 


how easily effects like equilibrium dissociation and 
jonization can be included in this solution, and they 
drew attention to the fact that a study of this type of 
flow can be very helpful in understanding the flow in a 
high-speed laminar boundary layer. 

In the present report, the Couette flow problem is 
extended to include a uniform externally applied mag- 
netic field normal to the walls. In the case of the sliding 
evlinders, the magnetic field will be radial (see Fig. 1). 
Unless stated otherwise, the magnetic field will be at- 
tached to the stationary wall. This is the case that is 
most interesting practically. 


GENERAL EQUATIONS 


Before writing down the equations that govern the 
Couette flow problem, it is helpful to write down some 
of the general equations and discuss briefly certain 
features of these equations. The steady continuum 
flow of a viscous, compressible, electrically conducting, 
electrically neutral (i.e., charge free) gas is considered. 
The rationalized m.k.s. system of electromagnetic units 
is used. 

From Maxwell's equations there is obtained 


curl E = (1) 
div B = 0 (2) 
j = curl (By) (3) 
Ohm’s law gives 
j= o(E + VxB) (4) 


where the conductivity o is assumed to be a scalar. 


The x B term is an effective electric field arising from 
the motion of the conducting fluid relative to the mag- 
netic field. 

The fluid mechanical equations will now be discussed. 
The continuity equation has its usual form 


div (pV) = (5) 
A convenient form of the momentum equation is 
div (pVV — P—T) = (Ga) 


The dyadic product VV in the transport term denotes 
the second-order tensor with the components 


(VV) a = (6b) 


P is the usual fluid mechanical stress tensor with the 
components 


Pu = —pbi + Tix (Gc) 


where 7, is the viscous part of the tensor (see, e.g., 
reference 4, pp. 333 to 338). T is the Maxwell stress 
tensor with the components 


T; 


= (BB, + eH 
(1/2) 


wu) te |E|?] (6d) 


HYPERSONIC COUETTE 


OF 


FLOW 603 
Div T is equal to the ponderomotive force per unit 
volume j x B acting on the fluid due to the interaction 
of the current and the magnetic field. Inserting the 
jx B term in Eq. (6a) leads to the more usual form of 
the momentum equation 


div (pVV — P) = jxB (Ge) 
The energy equation can be written in the following 


equivalent forms: 


div {pV[k + (V2/2)] — eV = 

(j?/0) + V-(jxB) =j-E (7a) 
div + (V22)] — 

eVtg+Ex(Bu)} =0 (7b) 


Here, h is the specific enthalpy, g is the heat flux vector, 
and ¢ is the viscous part of the fluid mechanical stress 
(6c) |. 
and div [E x (B yu)| represents the rate of transport 
(per unit area) of electromagnetic energy out of a fluid 
element. 

In general, the boundary conditions to be imposed 
In the Couette 


tensor {see Eq. ( Ex (B u) is the Poynting vector 


on £ and B are not easily determined. 


flow problem, £ will be set equal to zero—i.e., there will 
be no self-generated or externally imposed electric field: 
The physical situation that corresponds to E = 0is 
discussed in the next section. This means that no 
energy can be added to or removed from a fluid element 
by electromagnetic means and the energy equation 
reduces to its usual fluid mechanical form. The effect 
of the magnetic field is then confined to retarding the 
flow and increasing the dissipation of the energy of 
motion into heat. 


EQUATIONS FOR COUETTE FLOW WITH FE = 


It is clear from the simple geometry and type of 
motion that the velocities in the y and z directions will 
be zero for Couette flow—i.e., v = w = 0-—and that 
nothing can vary in the x and z direction —1.e., = 
0 0: = 0. The continuity equation is then identically 
satisfied and the transport terms in the momentum 
and energy equations are identically zero. From the 
z component of the momentum equation it follows, as 
expected, that B, = 0, the subscript denoting the com- 
ponent of B. Eq. (2) reduces to 

dB,/dy = 0 
so that 
(S) 


B, = constant = B, 


where By is the uniform externally imposed magnetic 
field. 


With E = 


0, Eqs. (3) and (4) become 


(9) 


the 
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Fic. 2. Temperature as a function of enthalpy and pressure. 


where j denotes the s (and only) component of the cur- 
rent density vector in the right-hand coordinate system 
of Fig. 1. In this equation and in the remainder of this 
report u is taken to be constant and equal to its vacuum 
value. 

The x and y components of the momentum equation 
reduce to 


(dr/dy) + (Bo/u) (dB,/dy) = 
(dr dy) — oBy*u 


(dp dy) + (d dy) (B,?/2u) = 
(dp/dy) — cB 


0 (10) 


0 (11) 


The first forms of these equations come from Eq. (6a) 

and contain the Maxwell stress tensor components, 

while the second forms come from Eq. (6e). The two 

forms are readily interchanged by using Eqs. (9). 7 is 

used here to denote the shear stress component 1,,, 

this being the only viscous component that remains. 
The energy equation reduces to 


(d/dy) (ru — q) = 0 (12) 


q denotes the y component of the heat flux vector, this 
being the only component that remains. 

The thermodynamics that will be needed here is the 
same as that described in reference 3 so that only a 
brief description will suffice. Thermodynamic equilib- 
rium is assumed throughout. The probable validity of 
this assumption in boundary-layer problems will be 
discussed later. The equilibrium thermodynamic prop- 
erties of air, the gas used in all the calculations in this 
report, are now available in a number of reports for 
fr: wide ranges of temperature and pressure (see, e.g., refer- 
q ‘ ences 5-9). The tables and/or Mollier charts in these 
3 references contain the thermal and caloric equations 
of state [i.e., f(p, p, T) = Oand h = h(p, T)| and thus 
define completely the thermodynamic state of the gas. 
Then a knowledge of any two thermodynamic variables 
yields all of the others. In this report, as in reference 
3, h rather than 7 turns out to be the natural variable, 
the other variable being p. Thus, all quantities that 
depend on the state of the gas will be considered to be 
functions of h and p. A plot of T(h, p) obtained from 
these references is given in Fig. 2. h/ho rather than 


h is used in this plot for convenience in the calculations 
where hf ~ is the free-stream enthalpy used in the caley 
lations. 
The usual Newtonian expression for 7 will be assume; 
i.e., 


7 = n(du/dy) (13 


where n denotes the viscosity coefficient. 
Likewise, the usual expression for g will be assumed 


q = —k(dT/dy) (14 


Here, k is an effective thermal conductivity that com. 
bines the usual thermal conductivity and a thermal 
conductivity that can be defined for the heat transfer 
in the y direction that arises because the dissociation 
energies are transported by the mass diffusion that js 
in turn, caused by the variations in the gas composition 
(see, e.g., reference 10). The sign is chosen so that 
heat flows from higher to lower temperatures. As 
noted above, it is convenient to write g in terms of | 
and p. Thus, Eq. (14) can be rewritten in the follow. 
ing form, using the fact that h = h(p, T): 


q = (n Pr) |—(dh dy) + (Oh Op)r (dp dy)| (1 


Here Pr = nC, k is an effective Prandtl Number and 


(Oh /Op)r is to be obtained from the equation of state | 
of the gas. In reference 3 (no magnetic field) the pres- | 
sure was strictly constant and the second term did not 
appear. However, here a new feature arises, in that 
the induced-magnetic-pressure variations induce fluid | 
pressure variations [see Eq. (11)], and both terms in | 
Eq. (15) would presumably have to be carried in 
general. 

The transport parameters o, y, and Pr will now be | 
briefly considered. For a more complete discussion | 
see Appendix (A). In general, these parameters are 
all functions of h and p. 7 and o are given as functions 
of h and p in Figs. 3 and 4, respectively. Pr is assumed 
to be constant in the formulas and calculations that 
follow. However, it should be emphasized that it is 
a simple matter to include a variable Prandtl Number 
in this problem if it should be desirable, and the pro- 
cedure will be indicated in the next section. 

The two Eqs. (9) together with Eqs. (10), (11), (12) 
(13), and (15) constitute a system of seven equations 
to be solved simultaneously for the seven unknowns 
B,, u, j, t, p, gq, and h. Since there are six first-order 
ordinary differential equations to be solved [Eqs. (4 
(omitting 7), (10), (11), (12), (13), and (15)], s1 
boundary conditions are required. It should be noteé 
that, if E were not zero, more equations and boundary 
conditions would be involved. 

Five of the boundary conditions are the usual ones 
imposed by fluid mechanical considerations. They are 


at y=0: u=0 

h=hy, or q = Qu 
at y=6: uw =U 

A= hk. 

= Po 
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The sixth boundary condition is to be imposed on B, 
and requires some discussion. First, consider the phys- 
ical situation that satisfies the condition K = 0. This 
can be seen by referring to the sliding cylinder configura- 
tion in Fig. 1. It is clear that the axial velocity « and 
the radial magnetic field By will lead to a circumfer- 
ential current density given by Eq. (9). If 6 is con- 
stant and flow conditions are identical around the 
circumference, the distribution of electric current being 
driven around the annulus by the effective electric 
field o( 1’ x B) is the same at all positioas around the 
circumference, and, thus, no electric field will be self- 
generated. If, in addition, no external electric field is 
imposed, the current is driven only by the effective 
electric field o( Vx B) and it is evident from Eq. (4) that 


this configuration satisfies the condition = 0. Now, 
note that this configuration is essentially a multilayer 
coil (solenoid) of infinite length, the magnitude and 
distribution of the current being known, of course, only 
after the problem is solved. It is well known that such 
a coil induces a magnetic field only in the axial direction 
and that the induced magnetic field is zero at the outer 
edge of the windings and is constant everywhere inside 
the coil. (See, for example, reference 18.) Thus, the 
boundary condition on B, is 


Since the current has the same sign throughout the 
annulus, B, will increase monotonically from zero at 
the moving cylinder to a maximum value at the sta- 
tionary cylinder. 

To completely specify the problem it is necessary to 
make some remarks about the properties of the walls. 
Since the moving wall corresponds to the non-ionized 
free stream, its conductivity will be assumed zero. It 
is not necessary to specify w for the moving wall al- 
though the vacuum value is indicated. If the moving 
wall were an electrical conductor, currents and a mag- 
netic field would be induced in it and it would be neces- 
sary to specify u. The solution is unaffected by the 
properties of the stationary wall (or body). Since the 
magnetic field is attached to the body no currents can 
If the body is ferromagnetic, » and 
0 but the solu- 


be induced in it. 
B, will change discontinuously at y = 
tion for y > 0 will be unaffected. 

At this point, a number of properties of the flow can 
be deduced with little effort. 

Changing the sign of By changes the sign of j and B, 
but does not essentially alter the solution. The reason 
is that the magnetic terms in the momentum Eqs. (10) 
and (11) enter as products or squares of these quanti- 


ties. Thus, j x By, which is the magnetic term in Eq. 
(10), always opposes the fluid motion, whatever the 
sign of By. Likewise, 7 x B,, which is the magnetic 
pressure term in Eq. (11), always alters the fluid pres- 
sure in the same way regardless of the sign of By. It 
should be pointed out, however, that there is a way in 
which changing the sign of By will alter the solution. 
This will occur if By changes sign in the s direction— 
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SONTC COVE 
i.e., circumferentially around the annulus—because 
then the current flow around the annulus changes sign 
and an electric field will be generated. Note that 
changing the sign of By in the x direction will not affect 
the solution. 

The first forms of Eqs. (10) and (11) can be integrated 
to give 


7 + (BoB, uw) = constant = 7, (16) 


b+ (17) 


2u) = constant = p, 
using the fact that B, = Oat y = 6. These equations 
show that both the shear stress 7 and the pressure p 
decrease monotonically from the free stream to the 
The well-known possibility of “holding a con- 
’ by magnetic pressure is 


wall. 
ducting gas away from a wall’ 
indicated here. 

The total current / and the maximum value of the 
induced magnetic field are, using Eqs. (9), 


J = jdy = -(1 dB, 
0 (Bz) 


ou dy (18) 
0 


If the quantities under the last integral are made non- 
dimensional using maximum values, there is obtained 


(0 Omar) (u U)d(y 6) (19) 
0 
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For a given set of boundary conditions, ¢/ ¢,,, and u U 
must be functions only of y/6 and, therefore, the value 
of the integral is independent of 6. This means that 
the magnitudes of the total current and the induced 
magnetic field are determined by the product Bod for 
given boundary conditions. It is clear that the B, can 
be made small or large compared with By by varying 6. 
The actual relative magnitudes of B, and By in a given 
problem when the boundary conditions By and 6 are 
prescribed will depend, of course, on the other factors 


in Eq. (19). 


From Eq. (19) it is seen that B, will be small com- 
pared with Bo if o,,,.; is small enough, or, better, if the 
magnetic Reynolds Number ¢,,,;46 is much less than 
unity. It should be pointed out that, even when B, is 
small compared with By, it is not correct to neglect B, 
compared with By in the general magnetic force term 
jxB=jx(B)+ B, =j7x Bo + jx B,. (One might 
be tempted to do this in order to obtain a simpler prob- 
lem which is a close approximation to the original 
problem.) One cannot neglect the B, term because 


jx By and j x B, are normal to each other and are there- 
fore involved in different components of the momentum 
equation [see Eqs. (10) and (11)]. This means that 
jx B, must be compared with other terms in its momen- 


tum equation (i.e., dp dy) rather than with 7 x By to 
determine whether or not it is negligible. Further- 
more, it should be noted that setting B, = 0 and retain- 


ing 7 x By would be an inconsistent procedure because 


j would then be zero [see Eqs. (9)]; indeed, if B, were 
strictly set equal to zero all magnetic effects would 
disappear. This inconsistency is not evident if only 
the second form of Eq. (10) is considered since B, is 
apparently not involved. However, the equivalent 
first form of Eq. (10) clearly shows the difficulty. 


Since it is necessary to retain 7 x By in the present 
problem, it is necessary to retain B, in Eq. (10). How- 
ever, it may be possible to neglect B, in Eq. (11) if the 
free-stream pressure is not too low, regardless of the 
relative magnitudes o1 B, and By. In this paper it will 
be assumed that B, can be neglected in Eq. (11) and 
that, therefore, the pressure is constant throughout 
the flow field. This assumption is justified later by 
substituting the calculated results into Eqs. (11) and 
(15) and showing that the pressure effects are actually 


negligible. 


The foregoing remarks about the induced magnetic 
field are easily generalized. If the induced field is 
small compared with the applied field, the induced field 
component parallel to the applied field can be set equal 
to zero to obtain a simpler problem which is a close 
approximation to the original problem. Whether or 
not the induced field components normal to the applied 
field can be neglected must be determined by consider- 
ations similar to those discussed above for the B, 
component. 
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SOLUTION OF THE COUETTE FLOW EQUATIONS Wp; 
p = CONSTANT 


With p constant, the system of equations and the 
boundary conditions may be summarized as follows: 


j = —(\/u) (dB,/dy) = (2 


(dr dy) + (Bo/u) (dB,/dy) = 
(dr dy) oBy*u = 0 (2) 


(d/dy) (ru — q) = O (29 
= n(du dy) (22 
g = —(n Pr) (dh dy) (24 
At y=0: u =0 
h =hy or Gg = qu 
At y=6: u =U 
h =h, 
B,=0 


Taking p to be constant greatly simplifies the task ¢/ 
solving the Couette flow equations. Aside from re 
ducing the order of the system of equations by elimi! 
nating the unknown /p, the unknown B, is now mw. 
coupled from « and h. In other words, after subst 
tuting Eqs. (23) and (24) into (21) and (22), the mo 
mentum and energy equations can be solved for u an( 
h without explicitly involving B,. B, can then be ob- 
tained from Eq. (20). The coupling of all of the equa. 
tions when ? is not constant is due to the fact that 7 an/ 
o are functions of the unknowns / and /p, and p involve 
B, explicitly. In the present simplified problem th 
momentum and energy equations are still coupled sinc 
n and o are functions of the unknown h. 

When the flow is incompressible, 7, ¢, C,, and Pr ar: 
all constant, and the solution of this system of equ: 
tions is trivial. Now, all of the equations are esset- 
tially uncoupled. For example, « can be obtained fron 
Eq. (21) after substituting for 7 from Eq. (23), then! 
and B, can be obtained from Eqs. (22) and (20), respec 
tively. The results are useful since they can be com 
pared with the corresponding boundary-layer solution 
The results are listed in Appendix (B) and some con: 
parisons with Rossow’s results are made. 

Since n(h, p..) and o(h, p..) are not given analytically 
it is obvious that an “‘exact’’ solution must be carrie’ 
out numerically. It should be noted that attempt 
have been made to obtain relatively simple analyti 
solutions, by making crude but simple analytical ap- 
proximations to n(h, p.) and o(h, p.), but without! 
success. Since the numerical solutions are easil! 
obtained, it does not appear to be fruitful to attemp! 
more involved analytical work. It might be mentione’ 
that all of the numerical work was carried out by the 
author using slide rule and graphical integration. 

The purpose of the present section is twofold—firs! 
to present the equations for the quantities of interes 
in the forms that were found most suitable for the 1 
merical work, and second, to deduce, where possible 
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properties of the solution by examination of these 


equations. 


Asin references 2 and 3, the appropriate independent 
variable turns out to uw instead of y. It will be con- 
venient to write all of the results in nondimensional 


form. 
The second form of Eq. (21) can be written 


(dr du) [n(du/dy)| — noBo?u = 0 
= noBy*u du 


Integrating from the wall (r = 7, « = 0) and solving 


for 7, 
| 
T = Te + re’) nou du 
0 


This equation can be written in the nondimensional 


form 
Nu 


r 24 
2 (25 
\ 1+2 non di (25) 


where the following notation has been used: 


7 = 


“=u ey. n = F0mar 
= maximum conductivity when there is zero 


magnetic field 


skin-friction coefficient 


C, = te/(1/2)p..U? 


= 


free-stream Reynolds Number 
based on the length 6 


pol 6 Ne 


Nu = Bod V omar Nx = & Hartmann Number 


Fimaz 1S used because it is easily obtained and is con- 
stant for a given set of boundary conditions.  »,,,,, iS 
obtained from Fig. 4, using homez, Which may be cal- 
culated from the equation given in Appendix (C). Ex- 
cept for By, the zero subscript will always denote zero 
magnetic field. 4 and @ are functions of “# through 
h(a) = (hh) (a%), which must be determined before 
the integration can be carried out. 

It is clear that 7 > 1 and increases monotonically 
from unity at the wall (4 = 0) to a maximum value 
7. atu = 1. This was noted in the previous section, 
using a different argument. It is seen that By and 6 
enter the equation as the product By. This is a general 
result that can be deduced from the forms of Eqs. (20) 
and (21). 

To obtain an equation for the enthalpy as a function 
of velocity—i.e., h(u)—the energy equation (22) is 
used. One integration can be carried out immediately 
to give 


Tl — = (26) 


Substituting for 7 and g, using Eqs. (23) and (24), this 
becomes 


nu(du dy) + (n/Pr) (dh/dy) = —qe 
and can be further written as 


[d(u?/2)/dy|] + (1/Pr) (dh/dy) = —(qw/n) 
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To indicate how easily a variable Pr can be included in 
the results, it will be assumed, for the moment, that 
Pr = Pr(h, p.). Integrating the above equation, 


u? 4 dh dy du (27) 
— = = 2 


At the upper wall this becomes 


h co U 
(U?/2) + f (dh/Pr) = —Qe f (du'r) (28) 
hy 0 
The recovery enthalpy h, is the enthalpy at the wall 
for zero heat transfer there—1.e., h, = h, for g, = 0. 
From Eq. (28), 


hr 
f (dh/Pr) = U?/2 (29) 
For constant Pr 
h, = (h,/h.) = 1 + (PrU?/2h.) = 
1+ [(y — 1)/2]PrM.2 (30) 


where 7 is the ratio of specific heats for the free stream, 
AM. is the free-stream Mach Number, and it has been 
assumed that the free stream is cool enough so perfect 
gas formulas apply there. It is seen that the recovery 
temperature is not altered by the magnetic field. It is 
expected that the same result holds for the boundary 
layer just as in the incompressible case.' Note that, 
if pressure variations were included, there would be an 
additional term in Eq. (30) that would depend on the 
pressure gradient and then the magnetic field would 
alter the recovery temperature. For the cases calcu- 
lated here this effect is negligible. 

If gu is eliminated between Eqs. (27) and (28) there 
obtains 


h hr 
f (dh/ Pr) -f (dh/Pr) X 
hw hw 


u U 
if (du (du — (w?/2) (31) 
0 0 
For constant Pr this can be written 
u 1 
h = hy, + (h, — hy) if (dit (di /7) | 
0 0 
— 1)/2] (32) 


For zero magnetic field, 7 = 1 and the ratio of integrals 
issimply 7. With a magnetic field the ratio of integrals 
is always greater than 7 except at the end points where 
it coincides with #7. This can be seen in the following 
way. Since 7 increases with # the integral in the nu- 
merator falls below #, the deviation from # increasing 
with #7. Thus, if this integral were plotted against a it 
would be convex upward. Then, the ratio of integrals, 
which coincides with @ at the end points, is also convex 
upward and hence must be greater than 7. This means 
that for the same value of # the enthalpy / with mag- 
netic field is greater than the enthalpy /) without mag- 
netic field—i.e., h(i) > ho(ai)—except at the end points 


ows: 
(2) 
(9) 
(23 
A 
| 
| 
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where they coincide. Since 7#(¥) and i%(¥) are generally 
different, it cannot be concluded that h(¥) > /o(¥). 
However, it is clear that the magnetic field will increase 
the temperature over part if not all of the flow field. 
This result is, of course, expected from the general phys- 
ics of the problem. 

For an insulated wall (zero heat transfer), h, = h, 
and Eq. (32) becomes 


h=1+ — 1)/2] PrM.201 — (32a) 


which is the result for zero magnetic field. This means 
that = h(i), not that h(¥) = ho(*). The expected 
increased temperatures show up through the fact that 
and are different. 

An expression for the skin friction is obtained from 
Eq. (23), which can be written 


= ndu/7 


Integrating and writing in nondimensional form, this 
becomes 


(C,R; = n(dit 
0 


1 
Then (C,R;/2) n(du 7) (33) 
0 


and the velocity distribution is given implicitly by 


u 
j= a(dit 7) / i(dii/7) (34) 
0 0 


Comparing the skin friction with and without a mag- 
netic field, 


1 1 
C,R; (C,R5)o = n(du nodu (35) 
0 0 


It is convenient to retain R; so that different boundary- 
layer thicknesses can be accounted for in any compari- 
son with and without magnetic field. In general 
n(t) > jo(a%), except for the heat insulated wall, in 
which case 4(7%) = no(7), this following from the previous 
results for /(7%). Then, for the insulated wall the mag- 
netic field decreases C,R;. At hypersonic speeds h 
will be large enough over most of the flow field so that 
n(h) varies slowly (see Fig. 3). Then 7 is expected to 
increase much more than 4 and it is expected that the 
magnetic field decreases C,R; in general. This is the 
result obtained for incompressible Couette and bound- 
ary-layer flow. 

The total drag (stress) on the body is the skin friction 
plus the force (stress) due to the interaction of the mov- 
ing conductor and the magnetic field. One way to ob- 
tain the total drag is to integrate the force per unit 


volume j x By across the flow field and add 7, to it. 
An equivalent procedure is to note that the total stress 
is the fluid mechanical plus the Maxwell stress [see 
Eq. (16)] and is equal to 7... 7. is obtained readily 
from the second form of Eq. (21) by integrating from 


0 to d—i.e., 
fo = Be ou dy 
0 
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If the total drag (stress) coefficient is defined by 


Ca = 72/(1/2)p..U? 


the result, in nondimensional form, is I 30 
Nu 2) nou du 
25: 
Eq. (23) can be written in the form : 
Integrating from 0 to 6, writing in ‘ational 2.0- 


form, and comparing with the total drag with no mag 
netic field gives 


1 / 1 
0 0 


It is clear that the C,R; is always increased by a mag. 
netic field even though C,R;is reduced. Incompressib. 


Couette and boundary-layer flows yield the same result a 

The wall heat-transfer rate g, can be obtained from| 

Eq. (28). Taking Pr constant, using Eq. (30), ané 
solving for g,, yields O5- 

Gow = (ho — Pr f (du (38 
0 

Defining the heat-transfer coefficient by ok 
( 

Cr = Qu/(hw — hr) pol 
Eq. (38) becomes ma. 5 
for the he 


C= | Pr (du 
0 


This equation can be written in the following two non- 
dimensional forms: 


1 405 
= apr f (du | (39 
0 
J 1 3.54 
= ad 2Pr [di (r T @ ) (40 
\ 0 
The ratio of the heat transfer with magnetic field to 305 


that without magnetic field is then 


1 2.59 
(CiRs)o = (di = 
0 


1 l 1 205 
n(du »/f ina f (di/7) = 
0 0 0 


154 
[CaR; 4 f (41 
0 
It is seen that (C,R;)/(C,R5)o is always greater tha Or 
(C,R;) and always less than (C,R;) 
The third expression, containing the three integrals 054 


shows that the magnetic field decreases both the nu 
merator and denominator, that they are roughly the! 
same magnitude, and that, hence, the magnetic fiel’| 0: 
does not have a large effect on C,,R;. C 
Eq. (39) shows the important fact that Reynold: | : 
analogy—i.e., the proportionality of heat transfer at | 
skin friction—no longer holds when there is a magnet T, = 1,2( 


ine 
: 
i 
i 
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field. Incompressible Couette and boundary-layer 
flows show the same result. 
The induced magnetic field is obtained by integrating 


Eq. (20). The result is 
1 
= [BoRay (C,R; 2) | f (nan du 7) (42) 


where Ry = is a magnetic Reynolds Number. 
The current density is 


= (BoRy ou (43) 
and, using Eq. (19), the total current is 
F 


The procedure for solving these equations numeri- 
cally will now be considered briefly. A more detailed 
discussion is given in Appendix (D). The general pro- 
cedure is to solve Eq. (25), which comes from the mo- 
mentum equation, and Eq. (32), which comes from the 
energy equation, simultaneously, using an iteration 
scheme. The shear stress distribution 7(7) and the 
enthalpy distribution /:(7) are thus obtained and are 
then used to calculate all other desired quantities from 
the equations given in this section. For the case of 
the heat insulated wall, Eq. (32) reduces to Eq. (32a), 
and no iteration is needed since /(#) can now be ob- 
tained without knowing 7(7). 


CALCULATED RESULTS FOR AIR 


Using the data given in Figs. 2, 5, and 4 and the equa- 
tions of the previous section, calculations have been 
carried out for the following cases: (1) heat insulated 
wall for values of the free-stream Mach Number 7, 
from 10 to 30 in steps of 5, and (2) one case of heat 
transfer (7, = 1,200°K.) for /7., = 20 and 50. Other 
pertinent values used in the calculations are 


Po 


atm, 
y = 1.40, 


T.. = 300°K.., 
Pr = 0.70 
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The free-stream pressure corresponds 
100,000 ft. altitude. 

For convenience in the discussion of the results, it wij i 
be assumed that comparisons with and without mag. 
netic field are made for the same value of 6—1.e., R, = 
(R;)o. Then, skin friction, heat transfer, and 
drag can be referred to explicitly, rather than CR | 
etc. The results will, however, always be presented jy 
the more general form with R; exhibited. In the mks 
system, B has the units (webers meter’). However 
in the presentation of the results, the more familiar’ 
gauss unit is used [1 (weber meter?) = 10‘ gauss]. 

No attempt will be made in this paper to determin, 
what constitutes practical magnetic field strength 
from the design point of view. However, for discussios| 
purposes, a rough classification of what will be referre 
to as moderate or reasonable magnetic field strengths 
will be given. Magnetic field strengths of the order o/ 
500 gauss or less will be referred to as reasonable sinc 
they can be produced with relative ease over a fairh 
large area without the use of extremely heavy equip. 


Note that a steady field strength of 10* gaus) 


roughly 


Fic. 11 
for the h 
ment. 1,200°K., 


may be taken as a practical maximum in the laboratory | 


Heat Insulated Wall 


Although this case is probably not very practical 
because of the high temperatures at the wall, it ij 
interesting to examine the effects of a magnetic field an‘ 
to compare the results for this case with the results fe 
the heat-transfer case which will be given later. 

Only the results for /,, = 30 are presented, the re| 
sults being typical for all the Mach Numbers con! 
sidered except for the magnetic field strengths te) 
quired. Results for 1/.. = 30 are presented in Fig 
5, 6, and 7. Fig. 5 shows that a tenfold decrease i: 
skin friction and a threefold increase in total drag ar 
obtained with the very reasonable magnetic fie!) 
strength By = 17.5 gauss when 6 = 1| meter, and thestil 
reasonable value By 175 gauss when 6 = 0.1 mete, py, 13 
It is clear from Fig. 6 and Eq. (26) that the rate of worl heat-trans: 
done ru on a fluid element and the heat flux g out of the ieee 
element are decreased near the stationary wall ani 
increased near the moving wall by the magnetic fieli| 
Fig. 7 shows that the net effect is an increased temper’ 
ature across the flow field and an increased heat trans) 
fer to the moving wall. 

The large effects caused by relatively weak magnet’ 
fields can, of course, be traced to the relatively lar) 
electrical conductivities associated with the high te 
peratures. The maximum conductivity for the ca 
presented is about 1,600 mhos meter. This should & 
compared with 10’ mhos/meter for copper and 10°) 
mhos meter for a good electrical insulator like gla 
For MV. 20, the maximum conductivity is abo! 
170 mhos meter and Byé = 40 gauss-meter product 
the same changes in skin friction and total drag as me! 
tioned above for 1J,. = 30. For J/. 10 the ma\' 
mum conductivity is about 0.06 mhos meter and cot _ 
parable magnetic effects are obtained for Bod = 2,0 a 
gauss-meter. According to the rough classification u%f? = 10~* 
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p = atm). 


10 | 
6 Ca RS 
7 (Cr Rs), 
6 
4 
3 
2 
3 ° (GAUSS-METER) 
199 
63.2 
65.0 
| 85.2 
23 3 
102 85.2 
65.0 
£26 a 
85.2 fe) 
yal 


612 JOURNAL OF THE AERO/SPACE SCIENCES—OCTOBER, 19658 


in this report, the magnetic field strength required for 
MJ... = 10 is not considered reasonable. 


Heat-Transfer Case 


Results for 1/.. = 30 are presented in Figs. 8, 9, and 
10. Fig. S shows that, for the reasonable magnetic 
field corresponding to the value of Bod = 25, a 400 per 
cent increase in total drag is accompanied by only a 
20 per cent increase in the heat transfer. Fig. 10 shows 
the pronounced heating effect of the magnetic field, 
the maximum temperature being increased from about 
4,400°K. for Bos = 0 to about 5,600°K. for Boyd = 
The small effect of the magnetic field on the wall heat 
transfer shows that the temperature reduction due to 
decreased viscous heating near the stationary wall and 
the temperature increase due to increased heating near 
the moving wall are about equal near the stationary 
wall. These changes in viscous heating also explain 
the shift of the maximum temperature position towards 
the moving wall. 

Results for 1/,. = 20 are presented in Figs. 11 
through 16. 

Consider, for the moment, the skin-friction curve in 
Fig. 11. The solid curve represents calculations and 
the dashed curve is drawn in to illustrate the probable 
continuation of this curve. The peculiar hysteresis 
character of the result can be traced to the fact that, 
when there is no magnetic field, enthalpies at all points 
in the flow field are low enough so that o varies rapidly 
with /# (see Fig. 4). A simplified problem based on a 
crude model of this flow that retains the essential o 
variation with / has been studied and found to exhibit 
the same hysteresis property. The reason this be- 
havior is not found for the case 1/.. = 30 is that, there, 
the enthalpies at most points of the flow field are high 
enough that o varies relatively slowly with / (see Fig. 
4). Furthermore, this behavior does not occur for the 
insulated wall for any Mach Number because the maxi- 
mum conductivity is fixed for all values of Bod and this 
essentially restricts the variation of o for most points 
in the flow field. If the skin-friction curve is divided 
into three branches by the two vertical tangents to the 
curve, it is clear that the branch contained between 
these tangents must represent unstable flows. This 
means that, if Bod is increased from zero, the skin- 
friction will follow the upper stable branch until the 
vertical tangent point at 237 gauss-meter is reached. 
Further increase of Bod will cause the skin friction to 
drop abruptly to the corresponding value on the lower 
stable branch. Note that this hysteresis effect makes it 
possible to obtain large effects with reasonable magnetic 
field strengths. 

In Fig. 11 the short portion of one stable branch of 
the heat-transfer curve shows the heat transfer decreas- 
ing with increasing Bod. This may also be seen in Fig. 
14, where a portion of one curve has been dashed for 
clarity. The trend of the heat-transfer curve indicates 
that it will probably drop below unity for large enough 
values of 


In Fig. 12, C,,, denotes the heat-transfer coefficiey 
at the moving wall. Except for the temperature boung.) 


ary conditions, which do not essentially change th, 
problem, it is easily deduced from symmetry that cI 
and C,, correspond to the skin friction and heat trang. 
fer at the stationary wall when the magnetic field isl 
attached to the moving wall—i.e., when the magneti:! 
field is fixed with respect to the free stream. Further 
more, C, and C;, give the results for the moving wall 
Such problems have been considered in reference |, 
Note in Fig. 15 that the induced field 5, will be small 
compared with the applied field if 6 corresponds to re: 
sonable boundary-layer thicknesses of the order of ; 
meter or less. However, the same conclusion canng 
be drawn for 1/7. = 30. 


BOUNDARY-LAYER CONSIDERATIONS 


In this brief section a few remarks will be made abou: 
the expected applicability of these Couette flow result 
to laminar boundary layers. | 

One of the primary differences between a boundary 
layer and Couette flow is its freedom to grow. Boun¢-| 
ary-layer thicknesses with and without magnetic field | 
will in general be different. Rossow! has shown that! 
the incompressible boundary layer grows when the 
magnetic field is applied. The present results indicat: 
that one effect of the magnetic field is to increase the 
temperatures in the boundary layer and thereby further 
increase its thickness. Estimates of the boundary 
layer growth due to the heating effect alone indicate that 
R; (R;)o will always be greater than C,R; (C,R5)o giver 
by the Couette flow calculations. It is expected, then, 
that for a boundary layer a magnetic field will actually 
decrease the heat transfer, the decrease being estimated | 
to be of the order of 15 to 30 per cent. Furthermore 
the skin-friction reductions are expected to be greater 
than indicated for Couette flow and the increases in 
total drag, while still present, will not be as pronouncet | 
as in Couette flow. 

It is expected that the hysteresis effect will be presen! 
in boundary layers and that the strengths of the mag- 
netic fields to produce the effects will be of the order 
indicated by the present results. 


DISCUSSION OF SOME ASSUMPTIONS 


One assumption in the calculations is that the pres 
sure is constant throughout the flow field. This as 
sumption can lead to errors in two ways: (1) by chang: 
ing the values of 7 and o and (2) by modifying the heat 
transfer as shown in Eq. (15). 

Consider the first source of error, and note from Figs 
3 and 4 that relatively large changes in pressure woul! 
be required to significantly alter 7 and o. Using the 
values of the induced field obtained from the calcula) 
tions, the maximum magnetic pressures 2 
were compared with the free-stream pressure Pp. = 
10-? atm. In all cases the ratio of the magnetic t? 
free-stream pressure was found to be of the order 
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10-2 or less and, thus, the changes in pressure are truly 
negligible. Consider, now, what will happen if p.. 
is reduced, say, to the value of 10-* atm. Because of 
the relatively small effects on » and oa, the calculated 
results would be about the same as obtained in this re- 
port. Now, the maximum magnetic pressure is of the 
order of P.. for some cases, and reductions in fluid pres- 
sure to vacuum near the wall are indicated. For such 
conditions the analysis of this paper would clearly not 
apply. From such considerations it is expected that 
the present analysis cannot be applied for values of 
much below atm. 

~ Consider now the second source of error. Using the 
above results for the pressure changes, the ratio of the 
second term to first term in Eq. (15) was calculated and 
found to be less than 10~° for all cases. Here, again, 
reductions in p. will increase the second term (it 
actually increases both factors) and violate the assump- 
tion of constant pressure. 

Another assumption that has been used is that @ is a 
scalar and is independent of the magnetic field strength. 
For this assumption to hold it is necessary that the 
Larmor or cyclotron frequency w, of the electron be 
much less than the collision frequency we of the elec- 
trons (see, for example, references 11 and 12). Fora 
magnetic field strength of B = 1 gauss the Larmor 
frequency is a, = 2.8 X 10° eps and w, increases pro- 
portionally with B. Based on a collision cross section 
of 1O- em.* [see Appendix (A) |, the collision frequency 
for the pressure of 10~° atm. is about 4 X 10° and varies 
proportionally with the pressure. It is clear that if it 
is required that w, we < 0.1 then B must be of the 
order of, or less than, 10 gauss. Field strengths of this 
magnitude can always be obtained by making 6 large 
enough. Examination of the results shows that for rea- 
sonable boundary-layer thicknesses of the order of 1 
meter or less the condition w, < we will not be satisfied. 
If the properties of o are not drastically changed for 
#, ~ we then the present results will be applicable to 
boundary layers. It should be remarked that when w,; 
is not much less than we, o has a tensor character, the 
circumferential current will be decreased as compared 
to the current for a scalar o, currents other than circum- 
ferential will occur, and the conductivity will be a func- 
tion of the magnetic field strength. 

The last assumption to be discussed is that of ther- 
modynamic equilibrium. This means that the reac- 
tion times are short compared with the time for an air 
particle to traverse the length of the body and to diffuse 
across the boundary layer. The effect of diffusion 
across the boundary layer will not be considered. If 


lion. (ionization time) > firaverse (traverse time) 


then the boundary layer will be essentially an electrical 
insulator and magnetic fields will have little or no effect. 
If 


tion. > taiss, (dissociation time) 


then the equilibrium case indicates the maximum pos- 
sible magnetogasdynamic effects. Note, however, 
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that if is much shorter than both and the 
ionization and conductivity in the boundary layer 
will be greater than for the case of thermodynamic 
equilibrium and the magnetic effects would be en- 
hanced. The various possibilities are discussed because 
little is known of ionization times in air, and the last 
possibility is at least suggested by Blackman.'* 


PROBABLE EFFECTS OF SEEDING 


The term “‘seeding”’ refers to the insertion into the 
boundary layer of small amounts of some foreign ma- 
terial that ionizes at appreciably Jower temperatures 
than air—e.g., potassium. In this way relatively large 
electrical conductivities in the boundary layer can be 
obtained at the lower Mach Numbers. 

Clearly, the probable effect of seeding is to shift the 
magnetogasdynamic effects shown in this paper to lower 
Mach Numbers and to show larger effects at the Mach 
Numbers given here. 


CONCLUDING REMARKS 


To illustrate some of the implications of the results 
obtained here, consider, for purposes of discussion, a 
body traveling at hypersonic speeds. Let the body be 
slender enough so that the drag on the body is primarily 
due to skin friction. Let the average surface temper- 
ature of the body be maintained at a value of the order 
of that used in the heat-transfer case in this discussion 
ie., = 1,200°K. 

The drag must show up as heat produced by viscous 
dissipation. Part of the heat is found in the body and 
part is found in the wake of the boundary layer. If 
the total drag (skin friction) of this body is increased 
e.g., by increasing the speed—the heat transfer in- 
creases in proportion. The increased dissipation asso- 
ciated with the increased drag shows up again as heat 
in both the body and the wake. 

Now, consider what happens when the magnetic field 
isturned on. The total drag will increase and the mag- 
netic field acts asa brake. The results show that large 
total drag increases are obtained with relatively weak 
magnetic fields. At the same time, large reductions in 
skin friction occur and the heat transfer is practically 
unaltered. The total drag is primarily magnetic drag— 
i.e., the force produced on the body via the magnetic 
field due to the interaction of the electrically conducting 
boundary layer with the magnetic field—because of the 
large reduction in skin friction. Since the heat transfer 
is unchanged, the heating associated with the increased 
drag must show up entirely in the boundary-layer wake. 
Note that the effect of increasing the total drag without 
increasing the heat transfer can be obtained by attach- 
ing a drag producing device (e.g., a parachute) to the 
rear of the body. 

The ratio of the total drag coefficient to the average 
heat-transfer coefficient is an important parameter in 
these considerations. For the slender body with no 
magnetic field this ratio is of the order of unity (Reyn- 
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olds’ analogy) and, from the present results, this ratio 
is much greater than unity for the magnetic case. 

From these considerations, it can be inferred that if 
the magnetic and nonmagnetic cases are compared on 
the basis of the same total drag, the heat transfer for the 
magnetic case will be much less than for the nonmag- 
netic case. 

It is interesting to note that a very blunt body at 
high speeds is similar to the magnetic case in that the 
ratio of total drag coefficient to average heat-transfer 
coefficient is large compared with unity. The drag of 
the blunt body is primarily pressure drag and the dis- 
sipative heating is mainly in the wake of the bow shock 
wave. 


APPENDIX (A)—TRANSPORT COEFFICIENTS 


Electrical Conductivity, o 


The electrical conductivity of air has been measured 
by Lamb and Lin.'‘ In order to obtain values of con- 
ductivity for a wide variety of pressures, o was calcu- 
lated using the following simple kinetic theory formula 
(see, e.g., reference 14): 


= (e2/S~/m,kT ) (N./N) (A-1) 


e = charge ofan electron 

m, = mass of an electron 

k = Boltzmann constant 

T = absolute temperature 

S = averaged collision cross section for scattering 
of electrons 


N, = number density of electrons—i.e., number per 
unit volume 

N = number density of particles other than elec- 
trons 


All numerical constants have been absorbed into the 
unknown S. Sis a weak function of T (see reference 
14) but will be taken to be constant in this paper. 
Examination of the results of reference 14 shows that 
the value of S = 10~" cm.? will give sufficiently accu- 
rate results for our purposes. NV, and N were obtained 
from tables of the equilibrium composition of air given 
in references 5 and 7. 

Data for NV, has been given for 3,000°K. and up. 
The calculated values of o were extrapolated to lower 
temperatures in a reasonable way on a plot of o versus 
T. The dashed portions of the curves in Fig. 4 indi- 
cate the extrapolated values, the wiggles in the curves 
arising from the wiggles in the 7 versus / curves of 
Fig. 2 


Viscosity Coefficient, n 


Exact calculations of » at high temperatures are not 
available. To the accuracy of existing knowledge, 
Greifinger’ indicates that the Sutherland formula is a 
reasonable approximation for the range of temperatures 
of interest here. Dommett'® shows somewhat poorer 
agreement with Sutherland’s formula. Sutherland’s 


JOURNAL OF THE AERO/SPACE SCIENCES- 


OCTOBER, 1958 


formula is considered to be sufficiently accurate for the | 


present study, and the following form has been used:"’ | 
§ 


= + 114)/(T + 114)] (T/T | 
Prandtl Number, Pr 


References 15 and 16 and references therein indicate 
that for the range of temperatures occurring in this | 
paper the Prandtl Number is nearly constant. The 
Prandtl Number has, therefore, been taken constant 


here. 


APPENDIX (B)—INCOMPRESSIBLE COUETTE FLow 


The results for incompressible Couette flow with the 
magnetic field fixed to the stationary wall are summa- 
rized below. | 


a@ = sinh (Nyj)/sinh Ny (B-1) | 
C,R;/2 = Nu/sinh Nu (B-2) 
T, = 1+ (PrU?/2C,T.) (B-3) 
f= T, + (7, — T.)5 
(PrU?/2C,T.) (sinh? (Nn ¥)/sinh? (B-4) 


C,R; = 1/Pr (B-5) 


CiR;/2 = (C,;R;/2) cosh Ny = 
Nu cosh Nu ‘sinh Nu (B-6) 


B, /‘BoRx [cosh Nu — cosh (Nu¥) | 
Nu sinh Ny (B-7) 


= [sinh (Nu¥)/sinh Ny| (B-8) 


These results exhibit most of the qualitative features of 
the incompressible boundary-layer problem solved by 
Rossow.' Note here that the recovery temperature is 
independent of the magnetic field. Eq. (B-5) shows 
that the heat transfer is also independent of the mag- 
netic field for fixed R;. If Rossow’s result for the 
growth of the boundary layer is used, it is found that a 
reduction in heat transfer, approximately equal to | 
that given by Rossow, occurs. 

Rossow also studied the effects of a magnetic field 
attached to the free stream. The Couette flow results 
for this case again show qualitative agreement with 
Rossow’s results. The Couette flow results for a mag- 
netic field attached to the moving wall are summarized 
below. 


= 1 — {sinh — ¥)]/sinh Ny} (B-9) 
C,R;/2 = Nu cosh Ny sinh Ny (B-10) 

= 1 + (PrU?/CpT..) X 
[((Ni cosh Nj, /sinh Ny) — (1/2)] (B-11) 


T=1+4 [T, — 1+ (PrU?/2C,T.)] (1 — 3) 
(PrU?/2C,T..){sinh? [Nu (1 — ¥)]/sinh? Ny} (B-12) | 


C,R; = 1/Pr (B-13) 


The heat transfer is not independent of the magnetic 
field for fixed R; since the recovery temperature is not 
independent of the magnetic field. 
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the AppENDIX (C)—-Maximum ENTHALPY FOR ZERO 
sed: MAGNETIC FIELD 
(A-2)  —s- This result is easily obtained in an obvious way from 
Eq. (32) with 7 = 1. The result is 
licate Nomaz Ny + [2(y 1)PrM.*] x 
this | — 1)/2)PrM.? + 1 — hy}? 
The 
stant. APPENDIX (D)—NUMERICAL PROCEDURE FOR SOLVING 
| EQUATIONS 
Consider the simple case of the heat insulated wall. 
- For given boundary conditions, 4(#) is given by Eq. 
ithe | (32a). n(%) and &@(#) are then obtained from Figs. 3 
ima. | and 4 and the integration in Eq. (25) can be carried out. 


For a given value of 2Ny,/C;R;, 7(a) is given by Eq. 
(25). C,R; 2 can then be obtained using Eq. (33), 
and .V, and, thus, Boé, are finally obtained. The 
fact that Bod is obtained from the calculations explains 
the odd values shown in the Figures. 

Consider now the general case with heat transfer. 
To start the iteration, assume 7(7) (e.g., the distribution 
for zero magnetic field, 7(#) = 1, might be used). 
Then, using /i(%) obtained from Eq. (32), 4(@) and 
a(#) are obtained from Figs. 3 and 4 and the integra- 
tion in Eq. (25) is carried out. Assuming a value for 
2Ny C,R;, a new 7(#) is obtained from Eq. (25) and is 
substituted back into Eq. (32). The iteration is re- 
peated until two successive distributions 7(7#) agree to 
within the desired accuracy. It was found that if 
2N,, C,R; was held constant throughout a series of 
iterations, the iteration scheme did not converge. 
However, if 2N,,/C,;R; was modified at each iteration 
so that the value of 7 for some fixed value of 7 was held 
constant, the iteration scheme converged rapidly, 
never requiring more than three iterations. Byd is ob- 
tained from the resulting value of 2N;,/C,R; after 
C,R; 2 is calculated using Eq. (33). 
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The Transonic 


JAMES 


SUMMARY 


The general three-dimensional theory for compressible flow 
through axial compressor blade rows! ? is applied to transonic 
cases. Because of the occurrence in the transonic range of acous- 
tic resonance, the theory is modified to include first-order vis- 
cous and heat-conduction effects. This approach is adequate 
to bring most transonic compressor problems within the realm 
of linear treatment. 

Numerical examples of the behavior of a transonic blade row 
are presented and discussed. The results indicate a smooth tran- 
sition from the subsonic to the transonic, or mixed-flow, regimes— 
i.c., the degeneracy encountered in two-dimensional transonic 
applications of the small-perturbation theory is completely ab- 
sent. At the same time, the results emphasize the failure of any 
strip thecry, or ‘‘blade-element,”’ approximation in the transonic 
range and thereby cast serious doubt upon the usefulness of tran- 
sonic cascade studies. 


SYMBOLS 

a = speed of sound 

A = B’*(1 — h) = similarity parameter 

B = number of blades in blade row 

c(r) = local blade chord 

Cath) = axial projection of local blade chord 

é; = specific heat at constant volume, undisturbed 
fluid 

bn = specific heat at constant pressure, undis- 
turbed fluid 

Cp = (p — po)/(1/2)pxU? 

Ce = section drag coefficient 

Con = total wave-drag coefficient per blade 

Ce = total torque coefficient per blade 

Cy = coefficients for determining K,,;, 

g(x) = profile shape function 

h = ry/rr = hubratio 

k = thermal conductivity 

= 2nr7/B = blade spacing 

Maz = axial Mach Number 

Mr = resultant Mach Number at tips 

My") = resultant tip Mach Number for resonance of 
mode 

m,k = summation indexes 

n = mB 

p = static pressure 

po = static pressure of undisturbed fluid 

Pr = Prandtl Number 
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perturbation velocity vector = (1, 7%) 

Fourier expansion coefficients for radial thick 
ness distribution 

radius of tips = radius of shroud 

Yn = radius of hub 


g(x, 7, 8) = 
QnB = 


rs = radius of sonic cylinder 

R = gasconstant = ¢, — Cy 

Re = Reynolds Number = v 

t = ft; = time 

= temperature 

i = temperature of undisturbed fluid 

U = axial speed 

Avy, = discontinuity in velocity normal to mean 
(helical) surface of blades 

“1, 91) = space-fixed cylindrical coordinates 

(x, 7,9) = rotor-fixed cylindrical coordinates 

8 = (1 — M,,?)'" 

5 = ratio of specific heats, c,/c, 

€ = dimensionless acoustic diffusivity, Eq. (14) 

B, mw’ = shear and bulk viscosities, undisturbed fluid 

v = kinematic viscosity, undisturbed fluid 

p = density 

Pin = density of undisturbed fluid 

o = wr/U = dimensionless radius 

or = wi,/U = dimensionless tip radius 

os = wrs/U = B/M,, = dimensionless radius of 
sonic cylinder 

r(r) = local thickness to chord ratio 

To = thickness to chord ratio at hub 


1, 4,4) = perturbation velocity potential, space-fixed 
system 

perturbation velocity potential, rotor-fixed 
system 


i 


w = angular velocity of blade row, rad. per unit | 


time 


INTRODUCTION 


— in the transonic compressor derives from 


the experimentally observed fact that relatively | 


efficient operation of axial compressors can be ob- 
tained throughout the transonic regime, and from the 
circumstance that maximum thrust for a given jet 
engine requires maximum mass flow through the en- 
gine, which in turn implies transonic Mach Numbers 
at the diffuser exit. Consequently, the transonic 
compressor blade row, with a subsonic resultant Mach 
Number at the hub and a near-sonic or supersonic re- 
sultant Mach Number at the tips (see Fig. 1), is of 
considerable interest to the designer. 

The transonic blade row, thus defined, has not pre- 
viously yielded to systematic investigation. The 
usual approximation, which assumes a two-dimensional 
flow field in each elemental compressor annulus, is 10 
longer applicable because of the complicated, three- 
dimensional nature of the interference pattern betwee! 
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Fic. 1. Blade row geometry and notation, showing cylindrical 
coordinates fixed in rotor. Helical arrow indicates incoming 
flow viewed from rotor. 


blades and because of the strong influence on one an- 

other of different radial elements of the same blade. 
' Furthermore, no theory that is linear in the usual 
sense of acoustics—i.e., linearized by assuming small 
disturbances—exists for steady, two-dimensional tran- 
sonic flows. As a result, a whole region of the flow field 
cannot be treated at all by a first-order two-dimensional 
approximation. 

Yet compressors do run in the transonic regime, and 
in fact do so with quite acceptable efficiencies. De- 
signers of these machines, in the absence of a satis- 
factory theory to guide them, having relied upon ex- 
periment, using cascade tests in many cases. In the 
transonic regime, however, the results of such experi- 
ments are difficult to interpret, not only because of the 
theoretical difficulties described above, but because of 
the impossibility of setting up experimentally a truly 
periodic transonic flow in two dimensions. Applica- 
tion of the cascade experiment to actual blade rows in 
the transonic case, therefore, is open to serious doubt; 
the results of the present investigation tend to em- 
phasize this fact. 


THE GENERAL METHOD 


The linearized (acoustic) theory contains the ap- 
proximation that disturbances propagate relative to 
the undisturbed fluid at a sound speed appropriate to 
the undisturbed condition. Consequently, if this 
theory is applied to sources and sinks moving exactly 
at the unperturbed sound speed, disturbances will, 
according to the theory, accumulate. The accumula- 
tion is particularly severe in steady two-dimensional 
cases because of the reinforcement of each source by 
its identical neighbors, and in these cases the theory 
cannot be applied. In two-dimensional steady flow, 
the accumulation may be relieved by going to higher 
order theories that take into account, for example, 
propagation at the local speed of sound relative to local 
fluid velocities. 

There are important physical cases, however, in 
which relief of the accumulation occurs even within 
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the acoustic approximation. Examples that are well 
known include unsteady two-dimensional transonic 
flows as well as the familiar wing-sweepback effect in 
steady three-dimensional flows. The present study of 
rotating machinery in transonic operation will provide 
a further example of three-dimensional relief. Al- 
though one source (or radial element of a blade) may 
be operating at a resultant Mach Number of unity, its 
immediate neighbors are not, and their disturbances 
will not reinforce those of the sonic source. Because 
of this, no degeneracy of the acoustic theory occurs and 
a first-order flow pattern in the usual sense exists. 


For these reasons, the method adopted by the author! 
has been to abandon the two-dimensional approxima- 
tion at the outset and to study, to first order in small 
disturbances, the complete three-dimensional flow 
pattern of an axial compressor blade row of arbitrary 
geometry. This general approach, and its application 
to nonlifting subsonic and supersonic blade rows, was 
the subject of an earlier paper.” In that paper, and in 
reference 1, the way in which the two-dimensional 
theory can be extracted rigorously from the three- 
dimensional was investigated. It was found that so 
long as the resultant tip Mach Number is not too close 
to unity (it may be either subsonic or supersonic), the 
classical two-dimensional cascade theories can be ob- 
tained by taking the limit as the blade number and the 
tip radius simultaneously approach infinity while the 
hub ratio approaches unity. Relaxation of this limit 
leads to three-dimensional corrections to the usual 
strip-theory approximation. 

As one must expect, when the resultant tip Mach 
Number approaches unity, the above limiting process 
is no longer possible, and the first-order flow field is 
inherently three-dimensional. Consequently, one must 
make use of the complete three-dimensional formulas 
to obtain, for example, the first-order pressure dis- 
tribution and wave energy losses of a transonic blade 
row. Moreover, the occurrence of acoustic resonance 
when the tip Mach Number exceeds one must be taken 
into account and is handled within the present theory 
by including the effects of viscous damping in the field. 


Since the number of blades is large, and this number 
determines the order of the Bessel functions that occur 
in the general solution, advantage can be taken of the 
asymptotic formulas for cylinder functions of large 
order.'~* Asa matter of fact, it is this procedure that 
makes possible the extraction of the strip theory 
i.e., two-dimensional theory in each flow annulus with 
the radius as a parameter—from the general three- 
dimensional theory in subsonic and supersonic cases.” 


THE THREE-DIMENSIONAL SOLUTION 


Assuming small-perturbations and inviscid potential 
flow outside of thin boundary layers, one can construct 
the general three-dimensional solution for flow through 
axial compressor blade rows by superposition of cylin- 
drical wave functions [see Eqs. (13) and (15) of refer- 
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ence 2]. Basie to this analysis is the familiar acoustic 
wave equation. 

The general solutions obtained in this way are ap- 
plicable to blade rows of arbitrary geometry (within 
limits consistent with the assumption of small per- 
turbations) over a wide range of Mach Numbers. 
They are degenerate, of course, if the axial Mach Num- 
ber approaches unity, but no such degeneracy occurs 


g(7, 6, 


. 


mB~k 


where c(7) is the local blade chord and ¢ = wr/U. 

It will be seen that the contribution of each particu- 
lar solution (or mode) of the wave equation to this 
general solution is proportional to Q,,g—4(%)/Anp—.- 
The coefficients Q,,,—,(x) depend only on the blade 
geometry,’ whereas the quantities A,,,—, are of the form 


V — (or"/os?) (2) 


The eigenvalues A,,,—, are the solutions of the trans- 
cendental equations 


R'np(Kmp—x) 0 (3) 


which guarantee vanishing radial velocities at the hub 
and at the shroud. The functions R,,,2(z) are normal- 
ized linear combinations of Bessel and Neumann 
functions occurring in each particular solution and are 
defined in reference 2. 

Because the number of blades in an axial compressor 
blade row is generally large, we find ourselves seeking 
solutions to Eqs. (3) in terms of cylinder functions of 
very large order. In that case, advantage can be taken 
of the asymptotic formulas developed by Watsor* or 
Olver,‘ and this procedure leads to numerical tech- 
niques that can be greatly simplified through the use 
of a correlating parameter involving the hub ratio and 
number of blades.':* In fact, when 7 is large, the eigen- 
values A,,. can be written 


Ky2/n? = 1 + (2C,/n?”); 


and in the limit of very large m the coefficients C, be- 


come functions of 4m”? alone, where 


A = B*3 (1 — h) 


“ Examples of the coefficients C,(A) are provided in Fig. 
3 of reference 2. 


(n = mB) 


Resonance 


A very important feature of the general solution (1) 
is that under certain conditions one or more of the 
Amp-~ can vanish. This occurrence is a kind of acoustic 
resonance and requires special consideration. In fact, 
it is a phenomenon peculiar to transonic operation, and, 
therefore, the applicability of the present theory to the 
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if the relative Mach Number anywhere in the floy 


annulus becomes unity. 

We are interested here only in the latter case, fo; 
which the sonic cylinder 
the undisturbed relative Mach Number is unity—g 
within or near the flow annulus. This can occur only 
when the axial Mach Number is subsonic, and then Ey 
(13) of reference 2 applies—viz., for nonlifting blade 
rows, 


imB w 


0 m=1k=1 'r 


Cs" U 


Nox 


transonic regime depends upon our success in treating 
resonant modes. 

The phenomenon of transonic resonance was noticed 
by Davidson® but was not treated in his paper. That 
it occurs only in the transonic (mixed-flow) regime 
follows from the fact that if \,,, = 0 


or/os = Ky/n (n = mB) (4 


It can be shown® that A,, > » for n > 0, and, conse- 
quently, resonance appears only when the radius of 
the sonic cylinder is less than that of the shroud (and 
is really important only when it is actually within the 
flow annulus). 


MODIFICATION OF THE THEORY 


Resonant modes have the property of remaining un- 
changed in form upstream and downstream of the 
matching plane. Once excited, they require no energy 
addition to sustain themselves. Since the compressor 
blading nevertheless puts energy into these modes, no 
finite steady-state amplitude is possible in the absence 
of a method for draining energy away from them. 
If, however, a realistic mechanism, such as viscosity, 
for extracting wave energy from the system is included 
within the analysis, a finite amplitude can be estab- 
lished. 


Viscous Damping 


The original assumptions that led to Eq. (1) do not 
exclude the possibility of taking into account the first- 
order effects of viscous damping and heat conduction 
in the field. External to thin boundary layers, 4 
potential solution still has meaning since no mechanism 
for producing first-order vorticity is present except 
shear at solid boundaries. Hence, we may put 


q "1, th) 


to first order, even though we include small but finite 
viscous damping in the theory. 

In analogy with the treatment of one-dimensional 
sound wave damping given by Stokes and improved 
by Kirchhoff,® we may begin by linearizing the Navier- 
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Stokes and energy equations to obtain 


dg/0 = —(1/p.)Vip + ([(4/3)u + X 
grad, div; g — veurl,*g (5) 


(1/p.) (Op Oh) — (¥/ px) (Op/Oh) = 

(k/e,) (VPT/p.T.) (6) 
where Cr, Rk, Po, and p., all refer to the 
undisturbed fluid. yu’ is the bulk viscosity (assumed 
zero for monatomic gases) and is included because of 
its importance in the damping of sound waves’ when 
rotational degrees of freedom are present. 

Eq. (6) simply states that to first order a fluid par- 
ticle can gain or lose thermal energy only through heat 
conduction. The treatment given by Stokes neglects 
heat conduction and in effect assumes the Prandtl Num- 
ber = uc, Rk to be infinite. Kirchhoff took heat con- 
duction into account in his analysis.® 

Eqs. (5) and (6) can be reduced in conjunction with 
the linearized continuity equation 


(Op + divig (7) 
and the equation of state 
p = pRT (S) 


to give a modified acoustic wave equation for the per- 
turbation potential in the form 


= (1, a*) — (4/3) + + 
— 1)/Pr]} (v/a?) (0/04)Vi2¢ (9) 


where a? = y(p./p.). In the derivation of Eq. (9), 
terms of order v? have been neglected in addition to 
second-order terms in any disturbance quantities.* 
Eqs. (5) through (9) are all referred to a space-fixed 
coordinate system. 

The corresponding pressure formula can be written 
in the form 


b' = p — Po = —Ppads, + [(4/3)u + (10) 


in which the usual acoustic formula can be recognized 
when = p’ = O. 


Consistency 


It would be well to pause here to consider the con- 
sistency of the present approach in neglecting higher- 
order inertia terms [for example, the g:Vq term in 
Eq. (5)] while keeping first-order viscous terms pro- 
portional to v. Simple dimensional analysis shows 
these terms to be in the ratio of a Reynolds Number 
based on the perturbation velocity and the appropriate 
wave length of the problem (in this case, L/m for the 
mth mode, where L is the spacing of the blades). 

However, the question of consistency actually re- 
quires deeper analysis than the simple process, implied 
above, of taking ratios, since the terms in question have 
opposite cumulative effects on the fluid motion. It 
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is well known that the higher order inertia terms in the 
equations of motion represent “‘steepening’’ terms that 
tend to create shock waves as a given plane disturb- 
ance propagates away from its source. Sooner or later 
these terms would always cause the violation of the 
assumption of small-perturbations in plane flow if it 
were not for the operation of viscous damping in pre- 
venting the steepening. 

The acoustic (first-order) theory is generally con- 
sidered adequate, therefore, for any small disturbance 
that damps out in a distance short compared with the 
distance required for steepening of the wave into a 
shock. The question of consistency, then, becomes 
one of deciding which mechanism—that of steepening 
or that of damping—first becomes predominant. 

In the case of one-dimensional sound waves, both the 
damping distance and the steepening distance can be 
calculated. It turns out that the acoustic theory is 
more applicable to high-frequency disturbances than 
to low-frequency ones, so long as lagging internal de- 
grees of freedom other than rotation are not too im- 
portant.’ 

In the present case of cylindrical geometry there is 
evidence that the wave forms, particularly in the case 
of resonant modes, will not steepen at all. Maslen and 
Moore* considered waves in an infinite cylindrical 
tube and showed that if second-order terms are in- 
cluded in the calculations the resulting strong waves 
do not steepen into shock waves as would be expected 
in other geometries. For this reason the author is en- 
couraged to believe that the inclusion of viscous mech- 
anisms within the linear theory is, in this case, justi- 
fied. 

The above arguments pertain only to the consistency 
of including viscous terms in a first-order calculation. 
This is not to say that higher order terms, viscous or 
inertial, might not become important; indeed, we shall 
see later that certain isolated nonlinear effects will be- 
come significant in some practical cases. 


The Solution in the Presence of Viscosity 


In coordinates that advance (at a speed Ll’) and ro- 
tate (at a speed w, see Fig. 1) with the blade row, the 
flow is steady, and in these coordinates (7, @, «) particu- 
lar solutions of Eq. (9) can be obtained in the form, 
with g(r, 0, x) = O(n, 41, fh), 
= An exp ind + (w/l’) X 

((in Myz?/B?) + Ane’ (11) 
where the R,(z) are the cylinder functions defined in 
reference 2, satisfying Eqs. (3). 

The solutions (11) are identical with the solutions of 
the ordinary wave equation’ that were used to build 
up the general solution (1), with the exception that the 
quantities \,,. are replaced by 


= — ime [(8/os*) + (11/052) + 3] — (Km?/or’) + (1/os?) + (3/82) + 1}, 


>0O) (12) 


* For a rather complete account of this procedure, particularly as it applies to plane sound waves, see Lighthill.? 
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Dox’ = (Koy o7B) + (e/2) [1 (Ko? ‘o7 8") | (13) 
where 


= (4/3) + (u'/u) + 
[(y — 1)/Pr}} (or/os*) (1/Re); Re = Urr/v (14) 


Eqs. (12) and (13) hold, provided terms of order 
e’ are neglected. Moreover, because of the linearity of 
Eq. (9), the general solution can again be obtained by 
appropriate superposition of the particular solutions, 
and the result is simply given by Eq. (1) with \,, every- 
where replaced by },,.’ for n = mB 2 0. 

Typical values of e for turbomachinery in general 
use are on the order of 10~° to 10~’, and, in all such 
cases, the difference between Ao, and Ay’, as well as 
the viscous term in Eq. (10), can be neglected entirely. 
However, in transonic operation, when it is possible for 
one or more of the functions X,,. to vanish for 7 > 0, the 
importance of the correction terms in Eq. (12) is 
clear; the viscous terms in fact establish the steady- 
state amplitude of the resonant modes. Indeed, exactly 
at resonance we have 


This result is in analogy with the result obtained by 
Stokes for the damping rate of one-dimensional acoustic 
waves. We note that the amplitude of resonant modes 
decreases with the 3/2 power of the blade number, and 
we may conclude from this that the present analysis is 
meaningful primarily for configurations having a large 
number of blades. 

Even when the number of blades is large, however, 
exact resonance (A,, = 0) of, say, the fundamental 
mode m = 1, k = 1 may imply sufficiently large dis- 
turbance magnitudes for nonlinear effects to become 
important. This problem remains to be investigated, 
but the analysis of reference 8 suggests that a successful 
method of iteration to obtain higher order terms asso- 
ciated with a given resonant mode could be developed. 
Even so, this would still leave the problem of treating 
the possibly important coupling between the resonant 
mode and the nonresonant ones. 


Pressure Coefficients and Forces on the Blades 


The limitation of the present solution [Eq. (1)] to 
the nonlifting, or thickness, problem actually implies 
that the blades of the rotor (or turbine wheel) are what 
Busemann has called ‘‘limp bodies.’’ That is, the 
mean camber line of each blade adjusts itself to the 
local interference pattern induced by its neighbors in 
such a way that there is no lift. A useful compressor 
blade row, of course, provides lift, but through the 
linearization of the problem, the lifting and thickness 
problems can be studied separately. It is felt that 


because of the appearance in the transonic regime of 
wave drag associated with thickness, the thickness 
problem is not only of interest in itself, but also pro- 
vides an indication of events in the more general case. 

As pointed out by Rott,’ no wave-energy loss (and 
hence no wave drag) can occur until the resultant tip 
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Mach Number is greater than unity. Once o7 os >]. 
however, some of the functions X,,, will be imaginary. 
and a net loss of wave energy by acoustic radiation will ' 
occur. In that case there will appear a drag on the 
blades associated entirely with the structural thickness 
provided. It is one of the purposes of the present 
theory to compute this drag force and the resulting 
torque that must be supplied to turn the blade row, 
If we define the pressure coefficient 


C,(r, 0, x) = p’/(1/2)p..U? = 
—[2(u + ov)/U] (16 


we may calculate the section drag coefficient per blade 
from the formula 


Ca (r) = D(r) (1/2) pa 
1 
C,(r, wx/U, x) [Av,(r, (17 
0 


where ¢,,;(7) is the axial projection of the local chord, _ 
and @ = wx/U defines the zeroth blade surface. The 
total wave drag and total torque coefficients can be | 
obtained by appropriate integration of C,(r) over the | 
span.! 


NUMERICAL EXAMPLES* 


Specification of the Geometry 


The theory described above has been used to compute 
numerical examples of the pressure distribution over 
the blades of a nonlifting axial compressor blade row 
operating in the transonic regime. A blade row con- | 
taining 64 blades and having a hub ratio of 0.9 was | 
assumed. (The parameter 1 defined earlier is then 
equal to 1.6.) 

In order to simplify the computation of the coefii- 
cients Q,,2-;(x), the axial projection of the chord, ¢,,, 
was taken to be constant over the span of each blade, 
and the blade profiles at each radial station were as- | 
sumed to be similar. 

These assumptions lead to the relation, consistent 
with the present linear theory, 


Av, (r, x) = g(x)r(r) (18) 


where 7(7) is the local thickness to chord ratio, and 
7g(x)/2U is the rate of change with axial position of 
the half-thickness about the mean (chordwise) camber 
line. (It is merely a matter of convenience to write : 
the specification of the profile shape in terms of the 
axial position instead of the position measured along 
the helical blades themselves.) 

Consistent with the geometrical configuration treated 
in reference 2, 7(7) was assumed to be of the form 


t(r) = roh(rr/r) V (hor? + 1)/(o2 + 1) (19) 


where 7) is the thickness ratio at the hub. This dis- 
tribution has the advantage of providing a radially 


* The numerical results reported here were largely obtained 
using the IBM 650 Magnetic Drum Calculator at the Cornell 
University Computing Center. 
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Fic. 2. Transonic pressure distribution at Wr = 0.98 for 


an axial compressor blade row with parabolic-are profiles and 
B = 64,h = 0.9,(A = 1.6), L = 2Caz, Maz = 0.6. The solid 
lines give the pressure coefficients over the chord at five dif- 
ferent radial stations ranging from hub to tip. The broken line 
gives the corresponding results from two-dimensional cascade 
theory applied at the tips. 


decreasing thickness as well as thickness ratio, typical 
of usual design. 

With this particular geometry the coefficients 
QOng—.(x) can be obtained from the set of numbers 
Onp-.* (see Table 2, reference 2) through the relations 


(xX) = hroV h2or? + 1 g(x) (20) 


once the profile shape is specified. For parabolic-are 
profiles, such as those used for the subsonic examples 
in reference 2, we have 


g(x) = 4U[1 — (2x/caz)] (21) 


Finally, a value of the solidity ratio, L/c,;, equal to 2.0 
completes the geometrical specification of the examples. 


Operating Conditions 


In order to span the transonic regime, six different 
operating conditions for the above blade row were 
assumed. In each of these the axial Mach Number 
M.: was taken equal to 0.6, and the resultant Mach 
Numbers at the tips, 


Mr = MV or? + 1 (22) 


were given the values 0.98, 1.00, 1.033, 1.050, 1.071, 
and 1.090 to define the six different transonic cases. 
The first two Mach Numbers belong to what we shall 
call the “‘lower transonic regime,’ the middle three to 
the ‘‘mixed-flow regime.’’ The final Mach Number, 
1.090, yields a slightly supersonic relative Mach Num- 
ber at the hub, but the sonic cylinder is still very near 
the flow annulus. 
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The Lower Transonic Regime— My, = 0.98 and 1.00 


The first two cases, with resultant tip Mach Num- 
bers of 0.98 and 1.00, are quite similar and essentially 
subsonic in character. In both cases all the parameters 
\,, are real, the pressure distributions are symmetrical 
(for symmetrical blades), and no wave drag occurs. 
The pressure distributions given by three-dimensional, 
linear theory in these two cases are presented in Figs. 
2 and 3 as functions of x c,, at several different radial 
stations. At each Mach Number, singularities at the 
leading and trailing edges, typical of linearized sub- 
sonic theories, appear. In both cases the singularities 
are integrable. 

In the first case, with a tip Mach Number of 0.98, 
it is possible to compare the results with the pressure 
distribution according to the two-dimensional Prandtl- 
Glauert theory. It will be seen that, although the form 
of the distributions is similar, differences between the 
two are of such a magnitude as to emphasize the failure 
at this point of the technique developed earlier’ * 
of computing three-dimensional ‘“‘corrections’’ to two- 
dimensional theory. Three-dimensional effects take 
on vastly increased importance when the tip Mach 
Number approaches unity. 

There is, of course, no linearized theory of two- 
dimensional flow that can be used for comparison with 
the three-dimensional results in Fig. 3 (\/J7 = 1.0). 
The very existence of these results provides a satisfy- 
ing verification of the earlier conjectures of Rott,’ 
and of the statements made earlier in this paper regard- 
ing three-dimensional relief of the usual transonic de- 
generacy of the linear theory. The form of the pres- 
sure distribution given in Fig. 3 is such as to indicate a 
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Fic. 3. Transonic pressure distribution at Wr = 1.00, Mar = 


0.6, for the same blade row parameters as in Fig. 2. 
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TON YES 
F 0.975 
Ce 5 950 
Fic. 4. Transonic pressure distribution at M7 = = 1.0338, 


Maz = 0.6, for the same blade row parameters as in Fig. 2. The 
values of « [Eq. (14)] used for this and subsequent calculations for 
Mr > 1 correspond to Re = 6.5 X 108, Pr = 0.7, y = 1.4,’ = 
0. At this Mach Number the sonic cylinder is approximately 
midway between the hub and the tip. The broken lines repre- 
sent the pressure coefficients according to Ackeret theory applied 
at r/rr equal to 0.975 and 1.0. 
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Fic. 5. Transonic pressure distribution at My, = 1.050, 


Maz = 0.6, for the same blade row parameter as in Fig. 2. The 
values of Re, Pr, y, and pw’ are the same as in Fig. 4. 


smooth transition from the subsonic to the transonic 
regimes. 


The Mixed-Flow Regime 


The results for the mixed-flow regime (Figs. 4 
through 6) verify still further the absence of transonic 
degeneracy. In these cases, however, the phenomenon 
of transonic acoustic resonance comes into play, and 
the amplitude of resonant mode must be established 
by taking into account viscous damping as outlined 
earlier. As the Mach Number approaches the value 
for resonance of any particular mode, viscous damping 
affects, and finally controls, the amplitude of this mode; 
it turns out, however, that the viscous effects are im- 
portant only in a very narrow range of Mach Number 
for each mode. 

The values of ¢ used for these calculations are equiva- 
lent to a Reynolds Number of 6.5 X 10° with Pr = 0.7 
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and y = 1.4, if = 0. They range from 2.29 X 
for Mp = 1.033 to 2.47 K 10-7 for My = 1.090. 

In spite of the smooth transition from the subsonj 
to the mixed-flow regimes, the pressure distribution; 
plotted in Figs. 4-6 for Mr = 1.033, 1.050, and 1.07; 
indicate that several interesting and unusual effects re. 
main to be explained in this range. As will be dis. 
cussed below, essentially two phenomena predominate 
namely, the expected transonic acoustic resonance 
or, at least, near-resonance*—-and strong interferencg 
between blades. At the tip Mach Number of 1.03) 
the phenomenon of near-resonance is particularly sig. 
nificant, but it also contributes to the unusual appear. 
ance of the pressure distributions at 7 = 1.050 and 
1.071. To illustrate the importance of the near. 
resonant modes in the latter case, the pressure distri. 


butions at r/rr = 0.90 and 1.0 are plotted again in Fig. 7} 


with the nearly resonant modes suppressed—4.e., with 
the coefficients of the nearly resonant modes put equal 
to zero. 

At Myr = 1.033 the pressure distribution bears some 
qualitative resemblance to the Ackeret theory applied 
at r/rr = 1.0 andr/ry = 0.975. Interference from the 
rear portion of the blade just in advance of the one in 
question appears near the leading edge, as it would in 
the two-dimensional theory, and there is a general 
compression over the forward part of the blade followed 

* A mode is here considered to be ‘‘nearly resonant’’ when its 
\,,.? function is reduced by a factor of 10 or greater compared with 
the X,,.? functions of the next higher and lower modes at a par- 
ticular Mach Number. Exactly at resonance (A, = 0), neither 
the real part nor the imaginary part of the function dominates 
and the only drag associated with an exactly resonant mode 
is that due to viscous losses. However, a mode can be nearly 
resonant and still have a dominantly subsonic (A,, real) or a 
dominantly supersonic (A,, imaginary) character. This will 
have an important effect on the build-up of wave drag. A given 


mode is on the ‘‘subsonic side’”’ of its resonant condition if its 


\,,, is real and on the “supersonic side”’ if its A, is imaginary. 


-5 


Fic. 6. Transonic pressure distribution at Mp = 1.071, 
Mar = 0.6, for the same blade row parameters as in Fig. 2 
The sonic cylinder is approximately at the hub in this case, and 
the values of Re, Pr, y, and yw’ are the same as in Fig. 4. The 
broken line represents the pressure coefficient according t? 
Ackeret theory applied at the tips. 
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TRANSONIC 


by an expansion rearward. The much greater magni- 
tude of the pressure coefficients, especially near the 
tips, than those of the Ackeret theory is partly ac- 
counted for by the approach to acoustic resonance of 
the mode m = 1,k = 1. 

At My = 1.0331, four modes are severely resonated 
namely, the m-k modes 1-1, 5-2, 6-2, and 10-3), but 
the convergence of the series is so rapid that the 1-1 
mode dominates the final results. The near-resonance 
of the 1-1 mode, which has a nearly constant radial 
amplitude, 1s particularly interesting because it forces 
the pressure distribution over the inboard half of the 
blading to appear largely supersonic in character, even 
though this portion of the blading is operating sub- 
sonically. This happens because the mode 1-1! is on 
the supersonic side of its exactly resonant condition. 
In addition, the resonance of fundamental modes such as 
1-1 or 2-1 has important implications with regard to 
the build-up of wave drag as the tip Mach Number in- 
creases above unity. This will be discussed below. 

To study further the effects of Reynolds Number in 
the region of resonance of the fundamental mode, nu- 
merical calculations have been made at J/7 = 1.033 
for Re = 6.5 X 10% and Re = 6.5 X 10° (besides 6.5 X 
10°). The results show that at this Mach Number the 
effects of viscosity in limiting the resonance of the mode 
1-1 are very small when the Reynolds Number is of 
order 10° or greater, but that at Re = 6.5 & 10% the 
amplitude of the near-resonant mode, and the magni- 
tude of the resulting wave drag, are reduced signif- 
icantly. These results indicate the narrowness, at 
typical Reynolds Numbers, of the viscous-affected 
range of ./7, since exact resonance of the mode 1-1 is 
at Wr = 1.0292 (see later discussion.) 

Closer to resonance, viscous effects become more and 
more important. The amplitude at resonance of 
nearly all the modes, as established by viscous damp- 
ing, is sufficiently small to keep these cases within the 
realm of the linear theory. Nevertheless, a check on 


Fic. 7. Transonic pressure coefficients at the hub and tips 
for Mr = 1.071, Mar = 0.6, with resonant modes suppressed. 
rhis Figure is constructed from Fig. 6 by subtracting the contri- 
butions of the nearly resonant modes—viz., m-k = 2-2, 3-3, 4-4, 
‘-6, and 10-8, 
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the amplitude assumed by mode 1-1 at Wr = 1.0292 
shows that if the Reynolds Number is on the order of 
10° the nonlinear effects discussed previously must be 
accounted for, whereas at Reynolds Numbers on the 
order of 104 the linear theory remains sufficient. Since 
practical Reynolds Numbers for axial flow compressors 
in aeronautics fall in the range of 10% to 10%, it is clear 
that further study may be required at a few isolated 
conditions of resonance. * 

The fact that the pressure distribution at J/; = 
1.033 is forced by the transonic resonance of the funda- 
mental mode to take on an essentially supersonic char- 
acter even inboard of the sonic cylinder is a very strik- 
ing illustration of the failure in this range of any two- 
dimensional theory based on the strip-theory approxi- 
mation. When resonant modes are excited the condi- 
tions at a given blade section can be strongly dependent 
on conditions at all the other blade sections. No strip- 
theory, or “‘blade-element,”’ approximation is then 
possible. 

At Mr = 1.050 (Fig. 5) modes 3-2 and 10-5 are 
nearly resonant, while modes 2-2 and 7-4 are also 
weakly excited. These resonant modes contribute 
to the tortured appearance of the pressure pattern; 
however, in contrast to the resonance of the funda- 
mental mode, they do not force an essentially super- 
sonic type of pressure distribution at sections operating 
subsonically, both because they are higher modes (hav- 
ing from two to four nodes over the span) and because 
all of them except mode 10-5 are on the subsonic side 
of their exactly resonant conditions. 

Very roughly, then, the pressure pattern at \/; = 
1.050 seems to be one that changes from a predomi- 
nantly supersonic character near the tips to a_pre- 
dominantly subsonic character near the hub, and the 
changeover occurs in the region of the sonic cylinder. 
Superposed on this dominant pattern is a smaller-scale 
structure, or waviness, due primarily to the near- 
resonance of the modes mentioned above. 

In addition to the main features just described, inter- 
ference between blades, particularly interference be- 
tween different radial stations of different blades, can 
be recognized at this Mach Number. This phenomenon 
is clearer, however, at /7 = 1.071 and will be discussed 
in detail for that case. 

When the tip Mach Number reaches 1.071 (Fig. 6) 
modes 3-3, 7-6 and 10-8 become strongly resonant. 
In addition, modes 2-2 and 4-4 are only slightly less 
resonant than the above modes and also contribute to 
the unusual appearance of the pressure distribution. 
To illustrate their effect, Fig. 7 was constructed from 
Fig. 6 by removing the contributions of all the nearly 
resonant modes. Comparison of Figs. 6 and 7 will show 


* On the other hand, Lighthill* states that the rate of acoustic 
attenuation in air can be from 2 to 100 times that given by 
Kirchhoff, due to the presence of H.O molecules or other non- 
classical dissipative effects. This means that in actual practice 
the effective Reynolds Number for the dissipative process we are 
considering here can be much smaller and that the present linear 
treatment may be adequate for nearly all practical cases. 
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that the phenomenon of resonance satisfactorily ac- 
counts for the small-scale waviness of the pressure dis- 
tribution over the inboard section. On the other hand, 
resonance will not explain the large recompression re- 
gion over the entire rearward span of the blade. 


This recompression cannot be explained on the basis 
of two-dimensional theory, either; in fact, in the tran- 
sonic regime the strip theory is generally misleading. 
The Ackeret theory applied at the tips indicates that 
a disturbance from the leading edge of the tip of a 
given blade cannot interfere with the trailing edge of 
the blade in advance of it (Fig. 6). The Ackeret theory 
does indicate, however, that the root sections of the two 
blades will interfere in this way. Strip theory there- 
fore, predicts a recompression only over part of the 
blade span.* 


In order to understand the existence of the com- 
pression over the whole span, it is convenient to make 
use of the three-dimensional theory of Ordway,’ in 
which the criteria for interference between points in a 
rotating cylindrical system are made available. It is 
possible, as the Ordway theory shows, for the root sec- 
tion near the leading edge of a given blade to interfere 
with the entire rearward span of the blade in advance 
of it. This is because the root section is operating 
nearly at the sonic cylinder and disturbances from this 
region propagate along very shallow and _ severely 
distorted Mach ‘‘cones.”’ All this is symptomatic of 
the transonic regime, and close study of the Ordway 
theory clearly shows the failure of the two-dimensional 
approximation in this regime. 


In sum, the pressure distributions in the transonic 
(mixed-flow) regime are strongly affected by the phe- 
nomenon of resonance and by patterns of interference 
explainable only from three-dimensional analysis. 
Furthermore, one is led to conclude from the Figures 
that the approximation of two-dimensional flow in each 
elemental compressor annulus breaks down, and that 
neither cascade theory nor cascade experiment is safely 
applicable in the transonic regime. 


On the other hand, although some nonlinear effects 
remain to be investigated, the present linear theory can 
be regarded as largely successful in treating the tran- 
sonic blade row. Indeed, it has isolated for further 
study such nonlinear effects as do appear important 
in the mixed-flow regime at some Reynolds Numbers, 
and has provided a complete coverage of the transonic 
regime at other Reynolds Numbers. Its applicability 

*On the other hand, two-dimensional theory does predict 
interference over the whole span from the rear portion of the blade 
just in advance of the one in question, causing the expansion re- 
gion near the leading edge that is also obtained in the three- 
dimensional theory. It will be noted, in fact, that the Ackeret 
theory applied at the tips agrees very well at this particular 
Mach Number with the magnitude and extent of the interference 
region over the forward part of the blade. This close agreement 
must be counted as coincidental, however, since it arises from the 
cancellation of opposite effects. For example, when the resonant 
modes are suppressed (Fig. 7), the agreement is no longer more 
than qualitative. 
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is improved (within limits) by the existence in at. 
mospheric air of nonclassical dissipative effects.’ 


The Case Mr = 1.090 


In this case (Fig. 8) the radius of the sonic cylinder 
is slightly less than the hub radius and the resultant 
Number at the hub is greater than unity. Although 
this condition is not included in our definition of 
“transonic’’ (see Introduction) it is of obvious interest 
not only in itself but also as some indication of the way 
in which the theory makes its transition from the tran. 
sonic to the supersonic regimes. 

At this Mach Number, several important modes, in- 
cluding 1-2, 1-3, 2-3, and 2+, are still subsonic jn 
form——i.e., their \,,,-,’s are still real. The resonances 
of these modes remain to be passed through. Also, 
modes 7-7 and 8-8 are mildly resonant and contribute 
the small-scale waviness that is still present at this 
Mach Number. 

Aside from the effects of resonance, however, the 
presence of three-dimensional interference patterns, es- 
pecially over the rearward portion of the blades, 
dominates the scene in this case. Once again, the 
Ackeret theory predicts qualitatively the nature of the 
interference pattern over the forward part of the blades, 
but only three-dimensional analysis can explain the re- 
compression region over the rearward part. 

It is clear, therefore, that transition, within the 
theory, from the transonic to the supersonic regimes is 
far from accomplished at M7 = 1.090. The sonic cyl- 
inder is still too near the flow annulus, and some im- 
portant modes have not yet changed form. On the 
other hand, modes 1-2, 1-3, 2-3, and 2-4 will all 
have changed from subsonic to supersonic character 
by the time the tip Mach Number reaches 1.16, so it 
appears that more or less complete transition to the 
supersonic regime can safely be expected for My ~ 
1.20. When this transition is accomplished, the strip 
theory can once again be expected to apply with rea- 
sonable accuracy. 
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Fic. 8. Transonic pressure distribution at Wr = 1.090, 


Max = 0.6, for the same blade row parameters as in Fig. 2. The | 
radius of the sonic cylinder is slightly less than the hub radius in 
this case. The values of Re, Pr, y, and yw’ are the same as it 
Fig. 4. The broken line represents the pressure coefficient ac- 
cording to Ackeret theory applied at the tips. 
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Fic. 9. Examples of two possible drag-rise curves for dif- 


ferent Reynolds Numbers. Plotted are conjectured values of 
the total wave-drag coefficient (per blade) versus the tip Mach 
Number for a blade row of geometry similar to that of Figs. 2 
through 8. The broken line indicates the values, independent 
of Re, according to two-dimensional, linearized (Ackeret) strip 
theory. 


It should be mentioned that for tip Mach Numbers 
sufficiently great for all or most of the X,,, to be imagi- 
nary (say 7 = 1.10) the convergence of the series be- 
comes markedly slower than in the case for lower Mach 
Numbers. It is then helpful to make use of the tech- 
nique of averaging partial sums in order to improve the 
numerical accuracy of the results. 


WAVE DRAG 


Any finite drag obtained using the pressure coeffi- 
cients computed from the present theory is identifiable 
as wave drag in the usual sense, with the exception that 
wave energy contained in resonant modes is partly lost 
through viscous dissipation, a process allowed for 
through the use of Eq. (9) in transonic, resonant cases. 
Nevertheless, in such transonic flows, no effort will be 
made to separate the two types of losses; the total 
drag computed will be referred to as ‘“‘wave drag.” 

The pressure distributions plotted in Figs. 4, 5, 6, and 
8 have been used to obtain the total wave drag and 
total torque coefficients (per blade) at the correspond- 
ing tip Mach Numbers. These results, 


Cow Ca(r)d (“) 


r 
Q Vite! (“) 


Mr Cow Co/ 
1.033 1.540 1.196 
1.050 0.544 0.396 
1.071 0.305 0.233 
1.090 0.312 0.2355 


provide an interesting indication on the behavior of 
the drag-rise curve as the tip Mach Number increases 
above unity. They do not, of course, provide adequate 
numerical information for the complete determination 
of the wave-drag build-up, because of the severe rate of 
change of pressure distribution with Mach Number in 
this range. 

However, the four pressure patterns given, indicating 
as they do the importance of transonic resonance, sug- 
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gest that the drag-rise curve also will be importantly 
affected by resonance, and hence by Reynolds Num- 
ber! When a fundamental mode is nearly resonant at 
a supersonic tip Mach Number just below its exactly 
resonant condition (A,, = 0), it will contribute very 
little to the drag (by definition of resonance). But 
immediately above resonance it will change from a pre- 
dominantly subsonic to a predominantly supersonic 
character and will still be greatly amplified. It will 
then immediately contribute strongly to the drag. 

These considerations and the results already obtained 
lead to the conjecture that two typical drag-rise curves 
at different Reynolds Numbers might be of the forms 
shown in Fig. 9. Successive bumps in these curves 
correspond to successive modes becoming resonant. 
Further numerical study of this problem weuld be of 
considerable interest, and such a study should provide 
information that could be checked experimentally. 
Indeed, if the conjectured ‘‘fine structure’ of the drag- 
rise curve does in fact exist, experimental study of its 
magnitude and significance would appear to be im- 
portant. 

In Fig. 9, the drag-rise curves computed according 
to two-dimensional strip theory are plotted for com- 
parison purposes. Although the section drag at the 
sonic cylinder tends to infinity at the sonic cylinder if 
Ackeret theory is applied there, the resulting singularity 
is integrable and the fofa/ drag computed in this way 
is finite. 


THE MaAcH NUMBER OF RESONANCE 


The foregoing numerical results illustrate the im- 
portance of transonic resonance, as regards its effect 
both on the pressure distribution and on the wave drag 
associated with thickness. Indeed, they lead one to 
the conclusion that for axial compressors the really sig- 
nificant Mach Number in the transonic range is not 
My = 1.0 but rather the resultant Mach Number for 
which the fundamental mode m = 1, k = 1 becomes 
resonant. (The resultant Mach Number for which 
other low modes, such as 2-1, 3-1, 1-2, 2-2, etc., are 
resonant are also important, of course. Their values 
can easily be obtained from the formulas given below 
and from Fig. 3 of reference 2.) 
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Fic. 10. Plot of the coefficient C, (4 ) (cf. Fig. 3, reference 2) 


for determining the Mach Number for resonance of the funda- 
mental mode. 
(24) and (25). 


The parameter A = B*: (1 — h). See Eqs. 
For m + 1, replace B by mB and A by Am?/s, 
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In general, the resonance condition is given by Eq. 
(4) with mn = mB. Using the formulas for the eigen- 
values A,,,. in terms of the coefficients C,, and the fact 
that os = 6/M,,, we obtain for resonance, 


(| — = 2C,(Am*™) (23) 


where /;\"“) is the resultant tip Mach Number for 
resonance of the mode m-—k. 

In particular, we have for resonance of the funda- 
mental mode 


= V1 + (24) 


Values of C\(A) are plotted in Fig. 10 for .1 = BY’ 
(1 — h) ranging from 0 to 10.* It will be seen that 
these results form an almost universal curve for B > 32. 
Using this figure and Eq. (24), the values of My"? can 
readily be found for any given configuration. 

At M,°~" the wave drag of an axial compressor 
should, according to the theory, increase significantly, 
and J/;"~! is therefore an important design param- 
eter. For tip Mach Numbers below 1/77" or 
the wave drag should be very small; hence these tip 
Mach Numbers, lying between 1.0 and J/7,‘'~!, are 
also of interest. They may readily be obtained from 


the relations 
Mp) = Vi + (25) 


and Fig. 10. 

Another important fact that can be deduced from 
Fig. 10 is that for 4 7 2.5, resonance of the funda- 
mental mode m = 1, k = 1 will not, as in Fig. 4, force 
a supersonic-type pressure pattern over the whole blade. 
This is because, as C,(.1) approaches its asymptotic 
value for large .1, the phase angles 6,,z-; in the defini- 
tion? of the functions R,,z, approach zero strongly (cf. 
Figs. 3 and 4 of reference 2) and then each function 
Riv|Kneailr/rr)| is very nearly zero near the hub. 
Thus, for hub ratios smaller than that used in the pres- 
ent set of examples, resonance of the fundamental 
mode may not be so disruptive of the pressure distribu- 
tion near the blade roots. 


CONCLUSIONS 


A linear theory has been developed that can be suc- 
cessfully applied to axial compressor blade rows oper- 


* The value of 1/7" ~" in the present example (4 = 1.6) is 


1.0292. 
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ating transonically. Although the theory is equally , 
applicable to lifting as well as nonlifting blade rows | 
detailed study of only the nonlifting, or thickness 
problem has been carried out to date. Numerical ey. 
amples have been worked out, and these examples 
verify the success of the linear theory in the transonic | 
regime. Even so, they also indicate the possibility | 
that for some configurations certain nonlinear effects 
might be important when resonance of the fundamenty] | 
mode occurs. 

The numerical results illustrate the importance oj 
transonic resonance and lead to introduction of the 
resultant Mach Number for resonance as an important 
design parameter. They also emphasize the failure 
of any strip theory or “‘blade-element’’ approximation 
for axial compressor blade rows operating in the 
transonic regime, thus leading to the conclusion that 
neither cascade theory nor cascade experiment can 
safely be applied to transonic blade rows. 
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Thermal Stress in Curved Beams 


BRUNO A. BOLEY* ann EUVAL S. BARREKETTE** 


Columbia University 


SUMMARY 


Thermoelastic stresses are investigated in curved beams of 
constant properties and constant cross section under temper- 
atures varving in the radial direction only. Good agreement is 
found between stresses as calculated from elasticity theory and 
from strength-of-materials theory which assumes that plane 


sections remain plane. A comparison with the corresponding 


straight-beam formula is also included. 


INTRODUCTION 


b hes STRESS ANALYSIS of curved beams is based, in 
strength-of-materials theory, on the assumption 
that plane sections remain plane; when the depth-to- 
radius ratio is small, it is permissible to further simplify 
the analysis by disregarding entirely the effect of 
that is, by employing the elementary for- 
A comparison between 


curvature 
mulas of straight-beam theory. 
the results of these approximate analyses with the cor- 
responding ones of elasticity theory may be found, 
for example, in reference 1, for a rectangular beam with 
a circular median line under a pure bending moment; 


it shows that the strength-of-materials approach is 
quite accurate for most beams, and it indicates the 
limits of applicability of the straight-beam equations. 
It is the purpose of the present work to perform the 
analogous comparison for curved beams under the 
action of nonuniform heating. The principal results 
are presented and discussed in the next section of this 
paper; the derivations of the pertinent equations are 
given in Sections 2 and 3. Additional details may be 
found in reference 2. 


(1) PRESENTATION AND DISCUSSION OF RESULTS 


Consider a curved beam of arbitrary constant cross 
section and with constant properties, the centerline 
of which is an are of a circle lying in one of the cen- 
troidal principal planes (Fig. 1), under a temperature 
T varying only in the radial direction. An approxi- 
mate formula for the stress a, in such a beam may be 
derived on the basis of the strength-of-materials as- 
sumption that plane sections remain plane, and is 
(Section 2): 


oak = —T + (1/A) f TdA + (1/A) TndA — ff | — )/mo(1 + + 
no 


+o) f T dn — + Bm) /m(1 + Bn)] (1/A) f 


where a is the coefficient of thermal expansion, 
Young’s modulus, A the cross-sectional area. All 
unspecified limits of integration must be chosen so as 
to cover this entire area, and 


B = h R, y=2 h, (1 ‘B) R) 1] 


The dimension 7) is a function of the cross-sectional 
geometry and is defined by the equation 


S ly + y)] 
A 


(la) 


(1b) 


Yo + y)] = 0 


and therefore specifies the location of the neutral axis 
for a beam in pure bending at constant temperature 
(explicit expressions for ro for several cross sections are 
given in reference 3); all other geometrical quantities 
are defined in Fig. 1. 

The formula for straight beams, analogous to Eq. 
(1), 


—T + (1/A) f TdA + (h'n/I) Tn dA 
(2) 


Ea 
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Tin) an) | dA (i) 


where the moment of inertia / is 
l=h? 


Eq. (2) should provide a reasonable approximation to 
Eq. (1) for beams of small depth-to-radius ratio. In 
fact, the first term of the series expansion of the right- 
hand side of Eq. (1) in powers of 6 is identical with the 
right-hand side of Eq. (2). The expansion in question 
is given in the second row of Table 1. Further remarks 
concerning the accuracy of Eq. (2) for curved beams 


(3) 


are given later. 

In order to estimate the accuracy of the stresses cal- 
culated from Eq. (1), they are compared with those 
obtained from an equation based on elasticity theory 
for the special case of a beam with a rectangular cross 
section of inner radius a and outer radius b. For such 
a cross section the geometrical parameters for Eq. (1) 
are 
ry = (b — a)/In (b/a);— = [1/In (b/a)] 

[(b + a)/2(b — a)] 
and the elasticity solution for the stress a, is (see Sec- 
tion 3): 


(4) 


3 
a 
| 
mples 
ffects | 
ental | 
ce of | 
the 
2 
: 
= 
= 
5 
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oy = —(A,/r?) + (A2/a?) [2 In (r/a) + 3] + (2A43/a?) — Ea [7 — (1/r?) f Tr ar] (5) 
a 
where THI 
A, = (Fa/N) X STF 
120 a)*In (b/a) [21n (b/a) — 1] + (b/a)? -1 f T rdr —4(b/a)*1n (b/a) f Tr (r/a) dr | 
a a ¥ 
A, = (Ea/N) [2(b/a)? In (b/a) — (b/a)? + 1] Tr dr — 2[(b a)? — 1] Tr In (r/a) dr ¢ oy 
a a H 
J b b 
A; = (Ea/N) —2(b/a)? [In (b/a) |? f Tr dr + [2(6/a)? In (ba) + (b/a)? — 1] f Tr In (r/a)dr ( : 
a a ) F- 
u 
N = 4(b/a)? [In (6/a)|? — [(b/a)? — 1/7 J 
3 
u 
The right-hand side of Eq. (5) may also be expanded in e 
powers of 8, as shown in the fourth row of Table 1. : 
The first term of this expansion is identical with the CENTER LINE 5 é 
right-hand side of Eq. (2), for the special case of a rec- 2 
tangular section. 
The third and fourth rows of Table | represent the a 4 r 
expansions in powers of 8 corresponding respectively 
to Eq. (1) (for the special case of a rectangular section) axis; See €0.[1b] 
and to Eq. (5). It may be noted that the first and ELA 
second terms of these expansions are identical; it is L/ os 
therefore reasonable to expect that the two formulas Pia | | 
will yield very similar numerical results, especially J) Fa 
for beams with a small value of 8.* To make a quanti- , / 
tative comparison, numerical results were obtained : ; 
hoes Eqs. (1), (2), and (5) for three special temper- Fic. 1. Geometry and sign convention for curved beam. 
ature distributions—namely, a constant (7) = 7)), a 
linear variation = 7\s/a), and a quadratic vari- 
is quite accurate; the approximate maximum stress | 
ation (7 = 7 For a constant temperature, all 
? es differs from the exact one by about —4.9, —0.7, and_ 
three equations yield the same result—namely, o, = 0.; 
: ; —0.44 per cent of the latter, respectively for b a = 
Q-—-and all may therefore be considered to be equally ‘ 
: 3, 2, and 1.2 (i.e., for 8 = 1, 0.667, and 0.182, respec- 
valid. Since these equations are linear in 7°, this re- : 
tively). The location of the maximum tensile (for 
sult implies that the addition of a constant term to any Mig od 
‘ > 0) stress, as calculated from the two theories, is 
given in Table 4. In this case the straight-beam for. 0. 
a gg — mula does not give zero stress, but a comparison be- 
For a linear temperature variation, the result o, = 0 2 =r 
. i tween its results and those of the other theories is 
again follows from Eq. (2), but not from Eqs. (1) and wed ; . 
: : ‘ difficult to make because only with Eq. (2) does the 
(5); here the straight-beam formula is clearly not P ; : 
: ; : addition of a linear term in the temperature not alter 0 
adequate. Eq. (1) is, however, quite accurate, as : : 
= : : the stresses. If the linear term is small compared to 
may be seen from Fig. 2. The approximate maximum 
: the quadratic one, the straight-beam formula is rea- 
stress (at ry = a) differs from the exact one by about . 
’ Ss sonably accurate. For example, in the present case. 
—2.4, —0.8, and —0.06 per cent of the latter for b/a = See ‘ , 
(T = the maximum stress given by the) 
3, 2, and 1.2, respectively. The location of the maxi- : ‘ - : 
; ; straight-beam formula differs from the exact one by 
mum tensile (for 7, > 0) stress, as calculated from the : 
about +13.7, +8.7, and +1.5 per cent of the latter, 
two theories, is given in Table 3. An example, cor- 
: is Ee for b/a = 3, 2, and 1.2, respectively. Note that the 
responding to Eq. (5), of the stress distribution across ; P a 
straight-beam theory gives equal stresses at = 
the selection is presented in Fig. 3. -0: 
: ‘ +h/2, while the other theories do not. 
The analogous comparison for a quadratic temper- 7 : 
It appears from the above discussion that Eq. (! 
ature distribution is shown in Fig. 4. Again Eq. (1) 5 ae ‘ 7 
Eomhaduaees is likely to be satisfactory for most practical cases. 
* The expansions in powers of 8, analogous to those of Table It is worth noting that even in the case of a rectangular 
1, for a beam in pure bending at a constant temperature, are beam Eq. (1) is much simpler to use, from the stan¢- -0.: 
presented in Table 2. The elasticity theory and the strength- Fic. 
point of numerical calculation, than Eq. (5); for ex 
of-materials formulas from which these expansions were derived f aie . 
are given in references 1 and 3, respectively. It may be noted ample, for the quadratic temperature distribution, 
that the first term of each expansion gives the appropriate Eq. (1) reduces to 
straight-beam stress, and that the first and second terms of the 
last two rows are identical. The results of the thermal case are 09/EaT, = (h? 3a’) [1/(1 + Bn)] | —2n*8 = 
thus entirely analogous to those of the nonthermal one. 3n? + n[(8/6) — (1/12m) — (1/8)] + 4)} (6) 


| 
| 
ips" 
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TABLE I 
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EQ 298 
THEORY SECTION a, a, as 
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BEAM ARBITRARY | 
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Series expansions for thermal stresses in curved beams. 
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Fic. 2. Stresses in a curved beam under a linear temperature -00 
distribution, T = 7\(z/a), 7, > 0. 
-0.10 
Fic.3. Example of thermal stress distribution, 7 = 7)(2/a). 
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while Eq. [5] remains essentially of the same form as 
previously given. 
(2) DERIVATION OF Eg. (1) 


A view of a deformed element of the beam of Fig. 1 
is shown in Fig. 5. The new length, after heating, of 
a fiber originally of length [(7) + y)dé@] is 


IG +yt+ + | 
0 


TABLE 3 
Location of Maximum Tensile Stress for T = 7\(s/a), T; > O 
(7»/a) for (7m/a) for 
b/a Eq. (1) Eq. (5) 
3 1.839 1.842 
2 1.448 1.449 
1.2 1.097 1.096 
TABLE 4 


Location of Maximum Tensile Stress for T = 72(s/a)?, T. > O 


(7m/a) for (7m/a) for 


b/a Eq. (1) Eq. (5) 
3 1.946 1.960 
2 1.483 1.484 
12 1.099 1.098 
0.4 
Oe 
EQ Ts -MAX. TENSILE 
0.2 
b 
a 
O 7 


-0.2 \ 


-0.4 


EQ. [5] \ \ 

\ 
-0.8 \ 
-1.0 


Fic. 4. Stresses in a curved beam under a quadratic temperature 
distribution, T = 72(z/a)?, T2 > 0. 
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provided that plane sections remain plane. 
is thus 


The strain 


= [1/(n + y)] ro’ — aTdy+ 
0 


0 


Noting, however, that 


aTdy<y (7a) 
0 


and introducing Hooke’s law, one obtains: 
[1 ‘(ro + y)] — + aT dy + 
0 

(Ad6/d8) + y) = (o,/E) + aT (8) 
The two unknowns 7%’ and (Ad@/d@) are to be found 
from the equilibrium conditions 

JS = oydA =0 (9) 
The expression for the stress is then 


og = —EaT + [Ea/(ro + y)] X 


(»/f saa) if TydA — 
f (» T ay) laa +f dy) (10) 


Eq. (1) results immediately from this formula with the 
introduction of the notation of Eqs. (la) and with 


y=h(n-—m), = 0 (10a) 


(3) DERIVATION OF Egs. (5) 


The governing equation of the present two-dimen- 
sional plane stress problem in terms of the Airy stress 
function ¢ is: 


Vid = -aEV*T (11) 


(Continued on page 643) 
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The Boundary-Layver Equation for Avxially 
Symmetric Flow Past a Body of Revolution— 
Motion of a Sphere 


DAVID MEKSYN* 


SUMMARY 


The equations of the boundary layer for an axially symmetric 
i solid of revolution are reduced to a form identical 


flow past ¢ 
This method is 


with that derived by the author for plane flow. 
applied to the case of a flow past a sphere, and the frictional force 
and the separation point are evaluated. 

The equation is integrated asymptotically by the method de- 
veloped for plane flow, where the partial differential equation is 


reduced to an ordinary equation. A new straightforward method 


of step-by-step integration of this equation is used. 


INTRODUCTION 


HE EQUATIONS of the boundary layer for plane 

flow were transformed by the author! to an equa- 
tion in which the boundary conditions are independent 
of the velocity at the edge of the boundary layer. The 
equation contains a gage invariant function which is of 
decisive importance in connection with the phenomenon 
of separation; it is also easier to integrate than the 
usual boundary-layer equations. 

A short derivation of this equation was recently given 
by Goertler.* 

In the present paper an identical equation is derived 
for axially symmetric flow past a body of revolution, by 
making use of Mangler’s* procedure which transforms 
the usual boundary-layer equations for axially sym- 
metric flow to those of a plane flow. 

The equation is applied to find the friction and the 
separation point for a flow past a sphere, measured by 


Fage.4 


) THE BOUNDARY-LAYER EQUATION 


Consider an axially symmetric flow past a solid of 
revolution. The coordinates x, y are taken along and 
normal to the meridian curve 6, respectively, u,v are the 
corresponding velocities, and r is the radius of revolu- 
tion (distance from the axis). 

The boundary-layer equations of motion and con- 
tinuity are 


= 7 Ux) + v(O7u oy") | 
= 0 


u(Ou Ox) + v(Ou/dy) 


(0/Ox)(ru) + (0/Oy)(7v) 
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where ly is the velocity at infinity (directed along the 
axis of rotation) and ll) is the velocity at the edge of 
the boundary layer. 

Changing the variables, let 


where / is an appropriate length 
Mangler’s transiormation® is obtained from Eqs. 
(1.2) by putting U; 
For an arbitrary function /(x, 


B= Ui(r 
Uy’ \ 


u 9= 


(1.2) 


y), we have 


(Of/Ox), = (r2/I2)U (Of + 
(Of /08),(0/dx) [(Ur/Dy] = (?2/2)U, 
(Of/da) + (Of/0B) 


(Of/Ov), = (rl 


whence the equation of continuity is transformed to 


4) denotes differentiation 


(Om Oa) + (07/08) 


where (1.4) 


and the subscript x in Eqs. (1 
with respect to x, when y is constant. 
Transforming similarly the equation of motion, we 


find 


u(Ou/Oa) + 3(0%/08) = 
(Uy? — In U,?/da) + v(07%/087) (1.5) 
Introducing a stream function (a, 8), such that 
= Oy/08, 7 = —Oy/da (1.6) 
and new variables 
we find 
oy = L(Of/dc) = 
dy = (2vU ya)’ *[(f/2a) + (OF 2a) (1.8) 
(a/2a)(Of, 
The equation of motion becomes? 
(0*f/00*) + = — (Of/de)*?] + 
2al (Of/de)(O°f/ Beda) — (4.9) 


where —a(d In U?/da) 


The boundary conditions are 


{ 
af 
und | 
if 
(9) 
+ 
| = 
he 
| % 
2 
3s 
} 
|| 


¢ 
¢= 


f = = 0, at 
Oa 


= 0, at o =0 (1.11) 
(2) FRICTIONAL FORCE 
The frictional stress F is equal to 
F = p(0u/dy),-0 (2.1) 


Now, 
du/dy = 
whence, after a few transformations, 


p(Ou, Oy) 


Sphere 


In the case of a sphere we put / = 7, where ro is the 
radius of the solid boundary and 7 is the distance from 


the axis; whence, after a few transformations, 


Cr = sin 0/a'/?) (a/R'”) (2.3) 


where 


a= (0? OC") R= 2roUo, a= sin? 


0 
In the application which follows, 


5.048 X sin 


R = 157,000, Cpr = eo (2.4) 


a’ 


(3) THe METHOD OF INTEGRATION 


The method of integration of Eq. (1.9) has been ex- 
plained elsewhere ;? however, to facilitate the reading of 


this paper, we give a very short account of it here. 
To integrate Eq. (1.9), assume that 


I(o, a) = (a,0"/n!) (3.1) 


where /(c, a) satisfies the boundary conditions (1.10) at 
Substituting Eq. (3.1) 


o = 0 for an ordinary point. 


dy = (1/3)(6a2)', dy = 
ds = (6/a)'*| [(@ — 28\a + 28aa’)/45] + (91d2/648a")} 


dy = G/ay [(62A — 75d*)/504] + (9ad’/28) — (ada’/3a) — (935d4, 15552a‘)! 
ds; = (3/a*){ [(—17 + GA)A2/280] — (Badd’/14) + (3A2aa’/10a) + 
where primes denote differentiation with respect to a. 
4+) EULER’S TRANSFORMATION 
(4) S(x) = [x/(1 + 4x)] = x — + 16x? — (4.1) 
The aim of Euler’s transformation is to eliminate a ; ; ; 
singularity which does not belong to the function itself for eg 1. ? The function S(x) has no singularity on the 
but is introduced by the method of expansion. The positive vere, the rp aera has, however, a singularity 
at x = 1/4, and it is divergent at x = 1. We could 


singularity can be eliminated by using a suitable ex- 
pansion. 
Consider as an example 
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(1.10) 


For the separation point, the additional condition is 
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in Eq. (1.9) we find the coefficients a, in terms of a. = q 
and its derivatives with respect to a. 


Substituting then Eq. (3.1), instead of f(o, a), in the | 


left-hand side of Eq. (1.9), and in the right-hand side 


we integrate it as a nonhomogeneous linear equation | 


in 0° //00*, and obtain 
= ¢(c, a) exp | —F(a, 


where Fe. a) 


a)do \ 


¢(o, a) being a slowly varying function. 
To find ¢(¢, a), put 


¢g(o, a) = (b,0"/n!) (3.3 
n=0 
Combining Eqs. (3.1), (3.2), and (3.3) and comparing 
on both sides the terms in o, we find the b,,’s in terms of 
a and its derivatives with respect to a. 
The unknown function a is found from the boundary 
condition (1.10), 


= f, ¢(o, a) exp — F(a, a)tde =) 


To evaluate the integral, we put 


F(o, a) (3.5) 


and invert Eq. (3.5) by expressing o in terms of 7; 
whence 


da 
f'(e,a) =1= e a) —dr_ (3.6)* 
0 dr 
where 


do 
dt m=0 
Integrating in gamma functions, we find finally the 
differential equation for a—namely, 


d,.V[(m + 1)/3] = 1 (3.7) 
m= 
where I is the gamma function. 
The expression (3.7) is divergent; it is summed by 
Euler’s transformation. Only six terms are used in 


the calculations, and they are given here for the sake of | 


convenience.°® 


d, = —(3d2/Sa?) 


obviously partly eliminate this singularity by a differ- 
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BOUNDARY-LAYER EQUATION 
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6 
Fic. 1. l) against 6, I = inviscid flow; II = measured flow 
(Fage). 


ent kind of expansion—for instance, by expanding in 
terms of 


y = + x) 


This is Euler’s transformation, 
The calculations for the case x = 1 can be conveni- 
ently arranged in the following table: 


1—4+ 16 — 64 
—3+ 12 — 48 
9 — 36 


—27 


where each term, starting from the second row, is an 
algebraic sum of two adjacent terms on the line above; 
whence 
S(1) = ] 3 9 24 
this expression is still divergent. Repeating the trans- 
formation, we obtain S(1) = 1/5. 

In the case of an asymptotic expansion (semicon- 
vergent), Euler applied the transformation in a different 
way. 

We consider a certain number of terms (in our case, 
six), apply the transformation, and retain only the 
convergent terms. The transformation is then re- 
peated on the remaining terms, and so on—this method 
will be used in our case. 

Since we make use of a limited number of terms only, 
all terms, including zeros, have to be retained. We can 
start the transformation from the first, second, etc., 
terms; the best result is obtained (with an equal num- 
ber of transformations) with a smallest last term and a 
best convergence. In any case, results corresponding 
to different combinations should differ only slightly. 


(5) PRELIMINARY RESULTS 


The above equations will now be applied to the case 
of flow past a sphere, measured by Fage. The main 
dimensions are as follows: d = 6 in., the diameter of 


the sphere; LU’; = 50 ft./sec., the velocity of the air; and 
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R = dly/v = 157,000, the Reynolds Number. Fage‘* 
measured the pressure and the friction. The minimum 
pressure is at 6 S 72.5°, and the separation point is at 
@ = 83° from the forward stagnation point. 

The functions and change rapidly between 
the maximum of l, and the separation point; and the 
values of U’,” can only be considered as approximate. 

In Figs. 1 and 2 are plotted the velocity 1 at the 
edge of the boundary layer and the curve of \ against 
6. The difference in the shape of A for a theoretical 
(inviscid) flow and a measured one should be noted. 

For a theoretical flow \ increases monotonically— it 
has a maximum for a real flow—and the separation 
point is a short distance downstream of the maximum. 


(6) INTEGRATION OF THE EQUATION— FRICTION 


The method of integration of the nonlinear ordinary 
equation (3.7) is step by step. In reference 5, the cor- 
responding values of a and a’ were found by trial, but 
since this is rather a tedious procedure, we are giving 
here a short straightforward method of step-by-step 
integration. 

Let the value of a be given at some point a = ay; we 
have to find a’, whence a is found for the next close 
value of a. The procedure is as follows. 

(a) We substitute the value of a in Eq. (3.7), leaving 
a’ indefinite; (b) the series being divergent, we apply 
Euler’s transformation and, from the expression ob- 
tained, find a’; and (c) substituting the value of a’ in 
the series from which it was obtained, we determine 
whether it is convergent—if only part of it is convergent 
we apply the transformation to the divergent terms to 
find a new value of a’—and so on. 

An example of these calculations is given below. 
Consider 


6 = 65°, A = —0.250, dd/da = 1.32, 
a = 0.378, a = 0.7353 
whence from Eq. (3.7) we find 
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Fic. 3. Cpagainst I, evaluated (Meksyn); II, evaluated 
(Tomotika ), ——— measured (Fage). 
f'(@, a) = 1 = 1.3430 — 0.38751 — 0.0413 + 


(0.2296 + 0.4198a’) + (0.4311 + 0.15070’) + 
(0.1188 + 0.0498a’) 


Applying Euler’s transfo. mation we obta n 


1 = 0.6715 + 0.2420 + 0.0689 + 
(0.0202 + 0.0262a’) + (0.0295 + 0.0572a’) + 
(0.0566 + 0.078la’) = 1.0887 + 0.1615a’ (6.1) 


— 0.549, and the last three terms become 


0.0058 — 0.0019 + 0.0157 


whence a’ = 


This is a moderately good result. One could apply the 
transformation to the last three terms; this, however, 
reduces considerably the magnitude of the terms of the 
equation, which makes it less reliable and does not 
improve the convergence. Then for @ = 70° and a = 
0.472, we find that 


a = 0.753 — 0.094 X 0.549 = 0.700 


The general procedure is as follows. At the front 
part of the sphere up toA ~ —0.4 (A = —0.5 at 6 = 0) 
we can disregard the terms in Eq. (1.9) in curly brack- 
ets, and take the values of a from Hartree® of the 
abridged equation;’ this has been done for @ up to @ = 
50”. 

For 6 = 55° we assume a = 0.840, which is slightly 
larger than the value a = 0.825 from Hartree’s solu- 
tion, since the curve of a against \ for the complete 
equation lies above the corresponding curve for Falkner 
and Skan’s’ equation. One has to bear in mind that 
the terms in Eq. (3.7) which determine a’ [see Eq. 
(6.1)] are small at the front part of the body, and, thus, 
are liable to considerable errors. It is therefore more 
practical to take, at the start, tentative values of a, 
until the terms in the equation are sufficiently large. 

The calculations were continued up to @ = 75° in 
intervals of five degrees; at 06 = 76.5° the convergence 
became less satisfactory, as it is the general case near 
the separation. The calculations were not extended 


beyond this angle. 

In Fig. 3 are plotted the coefficient of friction Cp 
[see Eq. (2.4)] against @ together with Tomotika’s 
calculations® up to @ = 55° by the momentum method 
(for somewhat different pressure distribution), and 
Fage’s* measured values from this angle to separation. 
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The separation point is evaluated by the method ey. 
plained in Section (7); the curve from 6 = 
drawn tentatively. 


76.5° ig 


(7) SEPARATION 


The separation point is found by a different set of 
equations than (3.7) and (3.8), since f(o, a) starts with 
the term oa’. 


The procedure was explained elsewhere;° here we give | 


only the final expressions, for the sake of convenience. 


The position of the separation point is found from the | 


equation 


d,V[(n + 1)/4] = 1 


n=0 


f'(@, a) (7.1) 


(7.2 


The remaining coefficients from dy) to dg inclusive are 
equal to zero; d; can be taken as equal to zero. 

The separation point is a short distance downstream 
of the maximum of \; we try, accordingly, few points 
until we obtain f’(~, a) = 1. 

As an example of the calculations, consider the point 


0.450, dd\/da = —7.00, a = 0.74 


6= 83°, A= 


Substituting in Eq. (7.1) we find 


f'(2, a) =0+4+ 1.4558 + 0 + 04+ 0 — 6.1227 + 


0 + 0 = 0.3639 + 0.3639 + (0.2738 + 0.1820 —- 
0.0181 — 0.2187 — 0.4024) 


Applying Euler’s transformation to the terms in paren- 
theses, we find 


0.1369 + 0.1139 + 0.0820 + 0.0409 — 0.0071 (7.3) 


whence = 1.0944. 

The transformation was repeated starting from 0.2738 
since this leads to a better convergence and a smaller 
last term than any other combination with the same 
number of terms. The last term in the Eq. (7.3) 
should be considered as equal to the algebraic sum of 
the two final terms; this is the procedure adopted by 
Euler, when two terms of very different magnitude 
follow each other. 

The corresponding result for 6 = 83.5° is f’(@) = 
1.0124, where the final expression is 


f’(e) = 0.3432 + 0.3432 + 0.1287 + 0.1072 + 


0.0739 + 0.0287 — 0.0125 = 1.0124 
The separation point is therefore at 6 = 83.5°+. 

It may be noted that this method of finding the sepa- 
ration point cannot be used with the theoretical pres- 
sure distribution since \ in such cases increases mono- 
tonically. 


(8S) Discussion 


The curve of friction obtained in the present calcula- 
tions differs somewhat from Tomotika’s results; this 
(Continued on page 664) 
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ERIK MOLLO-CHRISTENSEN* 
California Institute of Technology 


SUMMARY 


The general problem of utilization of test data in flutter anal- 
ysis is discussed. Methods for evaluating changes in flutter 
characteristics due to changes in structural and aerodynamic 
parameters are given. The precision of information obtained 
from tests is considered and the resulting accuracy of the flutter 


prediction discussed. 


SYMBOLS 


= streamwise coordinate 


v7 = spanwise coordinate 
o(x, ¥) = wing deflection or complex amplitude 


of wing deflection 


Poo = complex amplitude of impressed force 
(1/2)pl? = dynamic pressure 
M = Mach Number 


@ = frequency 


k = wh/l reduced frequency 
b = typical length dimension of wing 
1 = domain in (x, y) plane occupied by wing 
plan form 
w(x, y) = complex amplitude of downwash at wing 
surface 
D(k, = operator: w(x, y) = 
D(k, M) 2(x, y) 
C(x, ¥; 7) = flexibility function 
] = denotes a matrix 
= denotes a determinant 
m(x, ¥) = mass distribution of wing 
eulx, ¥) = pth natural mode amplitude distribution 
@(y) = uth natural frequency 


aerodynamic influence function: the 
lift coefficient per unit area at (x, y) 
due to a unit downwash at (£, 7) at 


reduced frequency k and Mach Num- 


P(x, 93 & 03k, M) = 


ber 
tn = generalized amplitude of the wth mode: 
n=(1 ff, Ax, y)dxdy 
M, = generalized mass of wth mode 
D = flutter determinant 
@ = subscript denoting test number 
[SHy | = diagonal matrix of integration weights 
Ou = generalized impressed force 
Quik, MW) = generalized aerodynamic coefficient 
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INTRODUCTION 


TREATMENT of the flutter problem as a part of 
design practice usually involves both flutter analy- 
sis and flutter tests. The latter are necessary because 
of the doubtful validity of the flutter 
cially at transonic and supersonic Mach Numbers. 
The discrepancy between results of analysis and the 
results are often difficult to explain, especially when the 
test observations are limited to flutter speed, Mach 
Number and frequency. 

It would be very desirable to establish a rational ap- 
proach to the problem of the discrepancy between 
theory and experiment as it occurs in practice, and in 
the process develop methods for extracting more useful 
information from flutter tests. Molineaux has sug- 
gested a method for finding the aerodynamic forces 
on an oscillating wing from flutter tests,' and Martuc- 
celli and the author have also outlined methods for 
this.* 

This paper discusses the general problem of obtain- 
ing and utilizing experimental data in flutter analysis. 
Since the practicable number of tests is limited, one 
first has to find the most important parameters of the 
flutter problem, and then use tests to determine those 
important parameters which are known from theory or 
previous tests with the least precision. 

The first step will be to write down the equations of 
motion. 


analysis, espe- 


(1) EQuaTIONS OF MOTION OF A FLEXIBLE WING IN 
AIR STREAM 
(1.1) Assumptions 
Assume the structure of the wing to be linear. Also, 


for convenience, assume the wing to be thin, so that 
the wing deflections are principally normal to a plane 
(x, vy) aligned with the remote airstream. It is then 
possible to use a flexibility influence function C(x, y; &, 
n), which is the deflection at a point (x, y) caused by a 
unit concentrated load at a point (&, 7). If one so 
wishes, one may consider the influence function to be 
modified according to reference 3, p. 109, so that 
rigid-body motions of the airplane are included. 

As far as the aerodynamic forces are concerned, as- 
sume the airplane to be in steady flight, which includes 
a steady maneuver, and consider the oscillations about 
the steady conditions sufficiently small to allow the 
equations of motion of the unsteady air motion to be 
linearized. These equations will be perturbation 
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equations upon the generally nonuniform steady flow, 
and will usually not have constant coefficients. Tran- 
sonic Mach Numbers are thus included, although the 
solution of the equations is not generally available, and 
one must, in order to obtain a flutter problem which 
can be solved, take recourse to brutal approximations. 

Neglecting such mathematical difficulties, there 
must exist an integrodifferential operator which gives 
the local lift at a point (£, ) for a deflection amplitude 
distribution z(r, s) for each Mach Number, M, and 
frequency of vibration w. If L(é, 7) is the local lift, it 
is expressible in the form 


L(é, n) = (1/2)pU? ff P(é, 8; k, M) X 
A 


D(k, M)z(r, s)dr ds (1.1) 


where 4 is the wing area, and k = wh/U the reduced 
frequency. For two-dimensional supersonic flow over 
a flat plate, for example, the operator in Eq. (1.1) is* 


L(g) = (1/2)pU"—4/V M2 — 1) (1/8) X 
[(0/0&’) + in) exp[—7kM?(é’ — r’) + 
0 


(M2 — 1)]Jo[RM(E’ — r’)/(M? 1)] X 
[((0/Or’) + ik] 2(r’)dr’ (1.2) 


where } is the semi-chord, and ~’ and r’ are obtained 
from &£ and r by division by 6. L(é) denotes the lift 
per unit span in this case. 


(1.2) Equations of Motion 


For a wing with a mass-distribution m(x, y), oscillat- 
ing harmonically in an airstream under the influence 
of impressed forces of complex amplitude F(x, y), the 
equation of motion can be written 


2(x, y) = ff. C(x, & [w’m(E, + 


F(é, n) + (1/2) pU? ff. P(é, 7, 5;k, M) X 


D(k, M)z(r, s)dr ds|dé dyn (1.3) 


For the solution of this equation in all but a few simple 
or simplified cases,®* ®one must take recourse to numeri- 
cal approximations. We shall consider two methods 
of approximation of the deflection amplitude distri- 
bution, z(x, y)—namely, (a) approximation by natural 
modes, and (b) approximation by discrete ordinates. 
This will in principle cover the forms of the equations 
resulting from any method of approximation. 


(1.3) Approximation by Natural Modes 


Let w,) and ¢,(x, y); u = 1, 2,... be, respectively, 
the sets of frequencies and amplitude distributions for 
free vibrations of the wing in a vacuum. Approximate 
z(x, y) by a sum of the first V natural modes: 

N 
¢, can be considered as a generalized coordinate of the 
wing motion. Substitution into Eq. (1.3) yields* 


1958 


N 


Where J/, is the generalized mass of the uth mode, 


M, = ff. m(é, n) [¢,(&, ) d&dn (1.6 


Q, is the impressed generalized force, 


If F(x, vy) y)dxdy (1.7 
A } 


and q’,,(k, MM) are the components of the generalized 
aerodynamic coefficients, given by 


= ff econ x 


P(x, y; & n) O(R, M) ¢,(&, n)dé ant dx dy (18 


Flutter is found as the condition for Eq. (1.5) to have 

a nontrivial solution in the absence of impressed forces. 
This requires the principal determinant of Eq. (1.5) to ~ 
vanish: 
(1.9 


| [(w? + (1 2)pU"q' M)]| 
where 


5 _ jl when p= vp 
“ when pp + 


The determinant is complex, and Eq. (1.9) can be solved 
to yield a pair of real characteristic values—for ex- | 
ample, flutter frequency w and air density p for a given 

Mach Number and atmospheric density and pressure | 
profile, which relates U’ to p and AV. 


(1.4) Approximation by Discrete Ordinates 


The integrals occurring in Eq. (1.3) can be approxi- 
mated by a sum over discrete values of the integrand. | 
The integral giving the aerodynamic coefficients in | 
Eq. (1.3) is observed to involve a differentiation 
D(Rk, M) 2(x, y). Although a derivative can be evalu- 
ated from discrete ordinates, the result is usually inac- 
curate. Partial integration with respect to the stream- 
wise coordinate may, if allowable, resolve this difficulty. | 
Alternatively, since the boundary conditions of the 
wing deflection are known, one may express 2(x, y) as 
a linear combination of slopes at discrete points 
(x,, as follows: 


N 
2(x, vy) ¥) (08/OK) + 
= 


o4(X, ¥) (02/Oy) eax, + + fo(y) (1.10) | 


y 


The functions \,(x, y) and o,(x, y) can be found for dis- 
crete values of x and y from experiment, for example. 
Introducing Eq. (1.10) into Eq. (1.3), one obtains a 
formulation of the problem in terms of slopes of the | 
wing surface. One may even go one step further and 
introduce second derivatives, which can be interpreted 
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as curvatures and twist, and, which can, incidentally, 
be measured using strain gages.” Using this, one may 
be able to obtain a good approximation to the elastic 
part of the problem and a reasonably sharp formulation 
of the boundary conditions of the aerodynamic problem 
using an approximation over a limited number of points 
on the wing. 

Whether the coordinates chosen are deflections, 
slopes, or curvatures and twists at a set of discrete 
points, the integral in Eq. (1.3) will appear as a sum, 
and Eq. (1.3) takes on the form 

\ 
+ F,) + 


N 
(1/2)pU? q,,(k, M)z,} (1.11) 
y=1 


where //, is the integration weight for the wth coordi- 
nate and gq,,(k, MM) are the modified aerodynamic 
coefficients. 2,;» = 1,2,...,/ N may stand for a 
deflection, slope, or curvature or twist, suitably ar- 
ranged in sequence. The condition for flutter is again 
given by the vanishing of the principal determinant of 
Eq. (1.11), which is the determinant of the matrix 


(—[1] + + 
(1/2)pU?[C] (1.12) 


Re(AD) 
Im(AD) 


Solving these equations for AU’, one obtains 
AU = 


or 


AU 


where the bar denotes the complex conjugate. It will 
be advantageous to perform these operations at con- 
stant Mach Number in an isothermal atmosphere, 
where Ll’, \/, w, and k are very simply related. Let 
p and T be the pressure and temperature, respectively. 
The velocity of sound, a, is constant during the change, 
and one has 


(1/2)pU? = y/2pM; U = Ma 
p = p/RT; w= kMa/b 
and the indicated derivatives are 
(OD = (OD/OM), (1/a) + 
{OD O[(1/2)pl2]! RT) + (OD/dw) (k/b) 
which is the change caused by a change in altitude, and 
(OD Ok)y = (OD/Ok) yy + 2[0D/O(w?) | /b?) 
which is the change caused by a change in frequency. 
Since the derivative of a determinant with respect to 


one of its elements is simply the cofactor of the element, 
these derivatives should not prove too difficult to evalu- 


(2) SENSITIVITY OF THE FLUTTER POINT TO CHANGES 
IN THE PARAMETERS OF THE PROBLEM 


The flutter determinant is a function of all the 
parameters of the problem. Call the determinant D. 
The flutter equation is 


D= D(k, p, M, U, m,..., Mn, 
Cn Cyn) = 0 (2.1) 


D is complex, while k, p, and .V/ are real. In addition, 
the relation between density, velocity, and Mach 
Number is given by the structure of the atmosphere or 
the characteristics of the wind tunnel, 


U = U(M, p) (2.2) 
The definition of reduced frequency k is 
k = wh (2.3) 


Buxton and Minhinnick have considered the rate of 
change of flutter speed with parameter changes,’ so 
have Baker> and Martuccelli and the author.’ Fol- 
lowing reference 2, one may observe that a change in a 
parameter P will change both the real and the imagi- 
nary parts of the flutter determinant. To bring the 
value of the flutter determinant back to zero, it will be 
necessary to change both velocity and reduced fre- 
quency, ./ and p being given by Eqs. (2.2) and (2.3). 

Thus, one has, as the conditions for maintaining 
flutter, 


AU Re(OD/OU), + Ak Re(OD/Ok)y + Re[(OD OP)AP| = 0 | (2.4) 
AU Re(OD/dU), + Ak Im(O0D + Im[(OD OP)AP| = Of 


—Re[(0D OP) AP |Im(OD/ok)U — Im[(OD OP) AP |Re(OD Ok) 
— Re(OD /Ok)y Im(OD OL’), 


to 


—Im {[(oD OP)AP] (OD dk) /Im} (OD (OD dk) } 


ate in practice. For the solution of a flutter problem 
on a high-speed computer, the elements of the flutter 
determinant will be stored in the memory, and there 
are standard methods for evaluating determinants, 
which only have to be modified by a program designed 
to skip selected terms. An alternative is, of course, 
to use numerical difference methods, involving repeated 
evaluation of the flutter determinant. 

The changes in flutter conditions which have been 
considered do only involve terms of the first order. 
Higher order terms may be considered, although this 
will increase complexity tremendously.’ Also, the 
method outlined above assumes that the flutter deter- 
minant is a continuous function of the problem param- 
eters. Although this holds true in most cases, it is 
perhaps not necessarily so, for example, at Mach Num- 
bers where a shock wave may jump from one place on 
the wing to another depending on boundary-layer 
characteristics or local separation. 
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In spite of this uncertainty, the method of analysis 
should have some practical value—specifically, it can 
be used to evaluate the effect of changes in wing prop- 
erties on flutter characteristics. It can also be used 
to find the precision of a flutter prediction, when the 
accuracy of the structural and aerodynamic informa- 
tion used can be estimated. Another application 
would be to enable one to identify a “‘sensitive’’ flutter 
point, where a small change in a wing parameter causes 
a large change in flutter speed. While the flutter pre- 
diction in this case is entirely unreliable, one has the 
comfort of knowing that it will not take much redesign 
to change the flutter speed. If, on the other hand, 
the flutter speed proves very insensitive to parameter 
changes, the flutter analysis will be very reliable, but 
it will take a major effort to budge a flutter point if it 
occurs within the flight envelope of the aircraft. This 
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state of affairs has been known in a less quantitative 
manner for some time, and can be used as an argument 
against radical innovations in a proven design. 

We shall here indicate how the method discussed jp 
this section can be utilized in connection with flutter 
tests. Before proceeding, we shall, however, outline 
a rough method of estimating changes in flutter speed 
due to parameter changes. This method will not work 
where one has a sensitive flutter point. 


(2.1) Rough Estimates of Flutter Point Changes Caused 
by Parameter Variations 


The method is based on iterative solution of the | 


flutter problem. The change in flutter speed is esti- 
mated by an additional step in the iteration, using the 
changed parameter. The method will thus only work 
when the iteration is strongly convergent. 

Eq. (1.3) can be solved by direction iteration" using 
the iteration formula 


2 (n) 
(x, y) = wi” ff Clx, é, n)m(&, 1) n)dédn + (5 l :) x 
A 


A A 


ww” [AM (x, y) + 1B (x, + (0 (C(x, vy) + iD (x, y)] (say) (2.6) 


with the normalization 
(Xo, vo) = (2.4) 


n refers to iteration number. From the normalization 


Eq. (2.7), one obtains: 


A Vo) 


B (xo, Vo) 
2)pl (x9, Vo) / Yo) 


D™ (x0, Yo) 
(2.8) 


(Xo, Yo), ‘B (x0, yo) | '2)pU?| (2.9) 


When this iteration process has been carried to con- 
vergence within a previously specified limit, which is 
accomplished in N, steps, say, one has determined 
flutter dynamic pressure and frequency. 

If one now wishes to find the effect of a small param- 
eter change, one may use y), and 
to carry out one more iteration, using a slightly changed 
value for the parameter in question, evaluating 
A™+)D(x, y), BU +) (x, y), ete. and substituting into 
Eq. (2.8). This will yield the new flutter dynamic 
pressure and frequency approximately, provided the 
convergence during the previous iteration was satis- 
factory. 


(3) EQUATION FOR A SET OF TESTS 


Because we shall consider the results of experiments 
as tools for determining the properties of the wing and 
the values of the aerodynamic coefficients, we shall have 
to use data from many tests to determine all the 
coefficients in the flight vibration problem. It will, 
therefore, be useful to write the equation for a set of 
tests as a single equation. To become accustomed to 


the thought, let us first consider the equation for a set 
of simple static tests. 

Consider the wing loaded with discrete loads IV’, at 
the points (x,, ¥,);# = 1, 2,..., Nand observe the de- 
flections z,, at the same points. Repeat the experi- 
ment NV times, and ‘denote the loadings in the ath 
experiment by W,, and the corresponding deflections 


by 
The equation of the first test is 
mer Ce... Cw | [Wa 
|= |]Cu Cr... Cw Wa 
(3.1) 
The equation for the N tests can be written 
Sul “pa “uN 
2NN_] Cy Cyn 
[Wu ... Wie... Win 
or, simply, 
[z] = [C] [W] (3.3) 


The matrix of flexibility influence coefficients |C] 
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can be found by post-multiplication of Eq. (3.2) by 
(Ww). 


iC] = fs) (3.4) 


Routine tests currently performed in industry avoid 
this last step by arranging the test such that the loading 
matrix {II’| is a diagonal matrix. Whether this yields 
optimum accuracy in [C| for a given accuracy of meas- 
urement of [2] and [I] will be discussed later. 

For other types of experiments, one may write down 
the equation for a set of N experiments in similar ways. 


(3.1) Equation for a Set of Ground Vibration Sets 


For a set of ground vibration tests, the equations of 
motion can be written 


= [C] [eye] + 
i[C] [Seas] + 
[C] [Free] 


(3.5) 


using discrete ordinates, where |d| is a matrix of damp- 
ing coeflicients, and a refers to the test number. Simi- 
larly, in terms of generalized coordinates, from Eq. 
(1.5), one may write: 


[SMA] [Spal tld’) [eax] 


where Q,. is the wth generalized impressed force in the 
ath test, and ay,) is the natural frequency of the uth 
mode, while w, (no parentheses around subscript) is the 
|d’| is a matrix of general- 


(3.6) 


frequency in the ath test. 
ized damping coefficients. 


(3.2) Equation for a Set of Flight Vibration Tests 


For sets of flight vibration tests the equations become 
more complicated unless one restricts the tests some- 
how. For example, one may consider tests carried out 
at constant Mach Number and reduced frequency, to 
obtain, using discrete ordinates, 


= [C] [Sm] + 


[C] [Fyal (3.7) 


where structural damping has been neglected. 
Formulating the flight vibration problem in terms of 
generalized coordinates [Eq. (1.5) | one obtains 


\f1 


(3.3) Use of the Equations for Sets of Tests 


The equations for a set of tests can be used to extract 
meaningful information from experiments beyond the 
usual observation that the test results do not agree 
with theory. Static tests have, of course, been used 
to obtain the flexibility influence coefficients for some 


time. One may, however, reason that the mass- 


distribution of a wing is known very accurately, so that 
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perhaps a ground vibration test would yield a more 
accurate matrix of influence coefficient than will a static 
test, since accelerometers can be used, rather than some 
reference plane which unavoidably will deflect. Other 
reasons for such a preference will be given later. 

One can also see how a set of model flutter tests or a 
flight vibration test can be used for the determination 
of the aerodynamic matrix [q,,(k, 17)|.. One may con- 
duct flight flutter tests and model flutter tests for such 
purposes, rather than interpreting test results using 
extrapolations based on generalizations from special 
It will first be necessary to consider methods 
Choose a 


~ases. 
for estimating the precision of test results. 
static test asa first example. 


(4) PRECISION OF A STATIC TEST 


Consider a set of static tests for the purpose of deter- 
mining the matrix of stiffness influence coefficients 
[C]. Assume that the loadings W,, are known ex- 
actly, and that the measured deflections z,, have the 
mean square errors o,*. It is, of course, open to ques- 
tion whether assuming that the errors are independent 
of the magnitude of the measured deflections is realistic, 
but the analysis will be simpler with this assumption. 

From Eq. (3.4) one obtains 


[C] = [2] (4.1) 


or, written out, 


\ 
s, = (W;,/|W]) (4.2) 


where J1’;, is the cofactor of and | is the deter- 
minant of the loading matrix [I1’]. 
The mean square deviation %c,;,; of C;,; is given by!! 


(o,2/|W|?) (4.3) 


One observes that, in order to minimize the error, one 


should make | IV’! large and >> H;,” small for all j, or, 
v=1 


stated differently, |I1’|~' should be small. The load- 
ings are limited to the elastic range, therefore one must 
seek other means to minimize the errors in |C|. The 
customary way is to make the loading matrix a diagonal 
matrix. Setting all the loads equal to unity, one obtains 


9 9 


= 


(4.4) 


as one should expect. 

If one chooses to express the equations of motion of 
the wing in terms or local slopes or curvatures and 
twists, the dispersion of errors in the appropriate in- 
fluence matrix will be given by equations similar to 
Eq. (4.4). 


(4.1) Errors in Mass-Distribution and Damping as 
Obtained From Ground Vibration Tests 


The purpose of ground vibration tests is usually to 
determine the natural modes of the aircraft in free 
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flight. Why this should be the purpose of the tests is 
not quite clear. It is impossible to suspend the aircraft 
such as to simulate free flight conditions, and the free 
vibration modes during the test are, therefore, differ- 
ent from the vibration modes in free flight. Also, one 
does not need to know the vibration modes in free flight 
to perform a flutter analysis. If a ground vibration 
test shall be useful beyond giving an indication as to 
which extent analysis is unreliable, it is necessary to 
consider it as a means for finding wing properties. If 
one is confident that one knows the mass distribution, 
one can compute the stiffness from the test results. 

On the other hand, if one is not confident about the 
numerical approximation involved in formulating the 
flutter problem in a finite number of terms, one may 
determine the constants in the equation of motion from 
tests, rather than analysis, and thus eliminate, within 
experimental accuracy, the remainder terms in the 
approximation as far as their effect on a set of specified 
vibration modes is concerned. In such tests, it will be 
necessary to measure all the forces on the wing or the 
airplane, including forces exerted by the supports. 

Consider the equation of a set of ground vibration 
tests, Eq. (3.5). Measure 2,,, # and the complex am- 
plitude of the impressed forces F,,,,, which includes forces 
exerted by the support. Assume that the matrix of 
flexibility influence coefficients [C] is known. 

Excite the structure as it is in one of its natural 
modes for each test by seeking a frequency where the 
amplitude is large and the phase angle between exciting 
force and motion is 90°. 

The deflection amplitudes will then be in phase, and 
we shall consider them real. From Eq. (3.5) one then 
obtains 


(4.5) 


for the real part of Eq. (3.5) and, for the imaginary 
part, 


[2] = — [wax X 


(4.6) 


The mass matrix [m] should, strictly speaking, be 
diagonal, and the magnitude of the off-diagonal terms 
as found from Eq. (4.5) gives an indication of the size 
of the mass-couplings which have been neglected. 

The precision of the mass matrix [m]| can be seen 
to depend upon the precision of [C|~! and of [z,.]7!. 
The errors in [C]~! are usually large, since the deter- 
minant of [C] is small for all practical wings. How- 
ever, the matrices [m] and [C]| only occur in the com- 
bination [C] [m] in the equations of motion, so it is 
more realistic to examine the errors in this matrix. 
From Eq. (4.5), 


To find the first-order error, assuming the frequencies 
to be exactly determined by the experiment, one differ- 
entiates Eq. (4.7) to obtain 
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d ((C] [Smy]) = (Sea? X 


(4.8) 


One observes that if [z,,|~' is very large, the errors 
will be large. Indeed, if two of the columns of [z,,) 
are very similar, [{z,,.|~! may be extremely large. Thus. 
if the structure has two modes of vibration closely simi- 


lar, only one of them should be included in Eq. (4.7), | 


If this is done, excessive errors will be avoided, since 
we have specified a test procedure which then yields an 
orthogonal matrix [z,,], which cannot be singular. 

The same remarks apply to the determination of 
[d]; also, it may be remarked that the accuracy in {d| 
depends directly upon the precision with which the im- 
pressed forces are measured. 

An alternative use of the data from a set of ground 
vibration tests is to determine the matrix of influence 
coefficients when the mass distribution is known. 


(4.2) Determination of the Structural Properties of the 
Wing From Ground Vibration Tests 


Perform the same set of ground vibration tests as 
specified in Section (3.1) above. When the mass dis- 
tribution [‘m,.] is known, the real part of Eq. (3.5) 
can be solved for [C], to yield 


Again, the errors in [C| behave as discussed in the 
previous section, but for the fact that they are multi- 
plied by [‘m,]~'. [‘m,] being a diagonal or almost 
diagonal matrix, this will not lead to difficulties. 

Solving Eq. (4.9) for [d| one obtains Eq. (4.6), as 
before, [C] now being given by Eq. (4.9). 

The reader is again reminded that [z,,| may be a 
slope or strain as well as a linear deflection with [C], 
[m | and [d] being interpreted as the appropriate general- 
ized wing properties. 


(4.9) 


(5) Use OF FLIGHT VIBRATION TESTS AND FLUTTER 
MopEL TESTS FOR DETERMINATION OF AERODYNAMIC 
FORCES 


The viewpoint adopted above can without modifica- 
tion be used for analysis of flight vibration tests or 
model flutter tests. For a set of N such tests, all per- 
formed at the same Mach Number and reduced fre- 
quency, Eq. (3.7) holds. It can be solved for [q,,(k, 
M) | to yield 


M)] = [C] [Pua] — [C] X 


(5.1) 


where [2,2], | pt) and [F,,] are 


measured. For a flutter test, F,, = 0, otherwise the 
equation remains the same. In terms of natural mode 
amplitudes, Eq. (1.9) can be solved similarly for the 

generalized aerodynamic coefficient q’,,(k, M) to yield 
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We proceed to discuss possible approximations for the 
aerodynamic coefficients, after which we shall return 
to the question of precision of aerodynamic coefficients 
and of flutter prediction. 


(5.1) Approximations to Aerodynamic Coefficients 


It will be inconvenient, if possible at all, to perform 
the entire set of flight vibration tests or flutter model 
tests at constant reduced frequency and Mach Number. 
It is therefore desirable to introduce approximations for 
the coefficients dependent upon reduced frequency and 
Mach Number. Realizing that, for high-speed air- 
craft with thin wings, the reduced flutter frequency is 
much smaller than unity, one may use an expansion 
in powers of k. Additionally, one may use some suit- 
able approximation for the variation of the aerodynamic 
coefficients with Mach Number. For supersonic 
speeds, it is possible to use powers of (J? — 1)'/?, 
while for transonic flow, powers of (J — 1) will be 
more suitable. For supersonic flow, the following ap- 
proximation is suggested : 
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dur (R, M) = M? 1) fa,, +, + 
(d,, + e,, (1/M?) +...4 (5.3) 
while, for transonic flow, 
Yur (R, M) = (ay, + dy (M 1) + 
+ d,, &) (M? — 1)? 4+... (5.4) 


where the a,,, b,,, etc., are constants. Better approxi- 
mations may, of course, suggest themselves. For the 
coefficients for some of the higher order terms in the 
expansions, theory may be entirely adequate, while, for 
the evaluation of terms which do not contain k, steady 
tests on deformed wings can be used. The remainder 
of the coefficients can then be found from flight vibra- 
tion or flutter model tests. The relative importance 
of the aerodynamic coefficients can be found using the 
expressions given in Section (2), and the unimportant 
coefficients can be guessed at, using available linearized 
aerodynamic theory. 

After all the aerodynamic coefficients have been ob- 
tained by experiment, guess, or theory, the flutter prob- 
lem can be solved, being careful of not exceeding the 
limits of validity of the approximations used for 
Gu(k, M). 


(5.2) Accuracy of the Aerodynamic Coefficients 


Differentiating Eq. (5.1), one obtains the first-order 
change in the q,,(k, 17) due to first-order changes in the 
terms on the right: 


a 


+ 


where it has been assumed that | (3 pl *) 


| [‘w,?, and are known exactly. To ensure that d [z,,‘]~! is 


small, the tests must be arranged such that [z,,.| is far from being singular—i.e., that the determinant of [z,,,] is large 


which can be done by exciting the natural flight modes. 
Appendix. 


Also, it is obvious that the errors in [z,,] and [‘m, | should be made small. 


An example of what otherwise might happen is given in the 


The error in [q,,(k, 7) | will 


be small when the dynamic pressures are high, which leads to the trivial conclusion that the most reliable flight 


flutter test is performed at flutter. 
hand side involving d[C]~!. 


The most bothersome term in Eq. (5.5) is, however, the first term on the right- 
As will be shown in the Appendix, this term is large for most structures, unless the prob- 
lem is formulated in terms of strains instead of deflections. 


There is still hope, however, since the aerodynamic coefficients only appear in the flutter equation in combination 


with the flexibility influence coefficients [C]. 
(C\[q,,(k, \f)]. The first-order error of this term is 


[qu(k, M)]} = — 


The term whose precision is important is, therefore, not [q,,(k, 47) | but 


([2yal — [C] [Fuel — [C} [Smy] |e (6.6) 


The important term in this relation is d[z,,]~', which, as has been mentioned before, can be minimized by testing 
at flight vibration modes, which will yield an almost orthogonal matrix [z,,]. 

At times, one may want to use the aerodynamic coefficients obtained from tests on one wing to carry out 
flutter analysis on another wing of the same plan form, but with different structural properties. 


Denoting 
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the matrix of flexibility influence coeflicients of the 
latter wing by [C’], the expression that occurs in the 


flutter equation is 


[C] — [C] |) X 


If [C’| and [C] are not very different, [C’] [C] will not 
differ much from the unit matrix, and the error will 
not be excessive. If, however, the wings are very dif- 
ferent, the error in [C’| [q,,| may become very large. 
For such cases, it will definitely be necessary to formu- 
late the problem in terms of strains instead of deflec- 
tions, since the strain influence matrices are usually 
further from being singular than are the deflection 
influence matrices. 


(6) AERODYNAMIC COEFFICIENTS OBTAINED FROM 
PRESSURE MEASUREMENT 


As an alternative to obtaining the aerodynamic co- 
efficients inferentially from deflection and or strain 
measurements, one may consider measuring pressure 
differences directly on oscillating wings. When doing 
so, the wing deflections must also be measured, other- 
wise the boundary conditions under which the meas- 
ured pressures occur would not be known. Also, the 
lift distribution is usually a more complicated function 
than the wing deflections and the wing strains, espe- 
cially in transonic flow. One must, therefore, use a 
large number of pressure pickups to obtain an adequate 
description of the pressure distributions. Also, in the 
author's experience with such measurements, shock 
waves seem to have a perverse affinity for the spaces 
between the pressure transducers. 

The wing strains, on the other hand, are given by 
weighted repeated integrals of the lift distribution, as 
are the deflections, and therefore much more smoothly 
distributed. 

Lift distribution measurements on oscillating flexible 
wings will therefore require much more effort than the 
aerodynamic admittance measurements described 
above, and should only be considered when one is 
mainly interested in throwing light on unsteady aero- 


AERO/SPACE 


SCIENCES—OCTOBER, 1958 
Further developments of traversing probe methods 
may modify this conclusion. 


(7) CONCLUSION 


Although most flight flutter testing will continue to 
be based on the generalization of trends found for par- 
ticular cases, it may be possible to use the viewpoint of 
this paper as a guide to explain the cases where accys. 
tomed methods fail. 

While no radical changes in test methods have been 
proposed, it is suggested that the purpose of the tests 
be viewed in another light, especially as far as ground 
vibration tests and flutter model tests are concerned, 
in the sense that tests and analysis be considered to- 
gether, and not as separate approaches to the flutter 
problem. The use, if any, of the methods outlined in 
this paper is best left to practising engineers. 

Evaluation of the methods which have been outlined 
can only be obtained by experiment, and such experi- 
ments are being planned, contingent upon the avail- 
ability of funds. 


APPENDIX—ILLUSTRATION OF THE EFFECT OF ERRORS 
IN A MATRIX UPON ERRORS IN ITS INVERSE 


As an example, consider the matrix of influence coef- 
ficients of a uniform cantilever beam under transverse 
loading, using three loading points—namely, at //3, 
21/3, and /—/ being the length of the beam. 

The matrix of deflection influence coefficients is, but 
for a constant, 


Cu Cre Ci3 2 5 Ss 
Cu Ce Cy} 16 28 
Cy Coo Co 8 28 54 
itsinverse, [A] = [C]~'is 
Ky Ky Ky 3.077 —1.768 0.462 
Ky Ky Kos | = | —1.768 1.692 —0.615 
Ky Ke 0.462 —0.615 0.269 
where A;,; is given by"® 
Ki = [C| (A-1) 


where C’ ,; is the cofactor of C;; in the transposed matrix 


dynamics, and not for routine practical measurements. [C;;|’. Differentiating Eq. (A-1) with respect to 
Table of (Cyr/Kij) (OK 
y=1 

| j=2 j=2 j=3 j=1 j=2 j=3 

| 
1 6.154 —6.154 —6.154 8.846 —8.846 8.846 —3.692 —3.697 3.692 
n= 1 2 —6.154 —16.413 —2.654 2.977 18.496 6.346 —6.266 —3.692 
3 —6.154 —2.654 — 1.582 —2.828 8.846 5.275 6.974 —1.192 —3.692 

= 
1 | 8.846 2.977 —2.82) —45.815 — 56.615 —45.949 7.431 12.361 17.231 
= 2 2 —8.846 18.496 8.846 — 56.615 —27.077 — 56.615 17.231 36.028 —17.238i 
| 3 8.846 6.346 —5.275 —45.949 — 56.615 — 52.044 5. 564 13.731 17.231 
1 —3.692 6.974 7.431 17.231 5.564 —3.739 —8.669 — 14.538 
w=3 |2 —3.692 —6. 266 —1.192 12.361 36.028 13.731 —8.669 —27.944 — 15.538 
3 3.692 —3. 662 —3.692 17.231 —17.231 17.231 — 14.538 —14.538 14.538 
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where C’,,,, denotes the cofactor of the (4, v) term in the 
cofacter of in the transpose of [C]. 

Eq. (A-2) gives the first-order relative change in 
[K,,|, (AK;,) Kj, due to a relative change in C,, of 
(AC,,) C,.. Table 1 shows these first-order changes. 

The numbers in Table | can be thought of at the 
matrix of percentage errors in [A | for a | per cent error 
in one of the terms in [C]. The reason why the errors 
are so large is that the determinant of [C] is small, the 
columns being almost simple multiples. Physically 
speaking, the deflection curves under concentrated 
loads are not really very different. The situation is a 
little better for most airplane wings, since their stiffness 
tapers off toward the tips. 

The inverse of the numbers in Table | are the first- 
order relative sensitivity of the terms in [C] to relative 
changes in [A]. These sensitivities are small. One 
must conclude that a given matrix of flexibility in- 
fluence coefficients does not specify the actual structure 
very precisely, while [A] seems to do so. 

When trying to evaluate [C] theoretically, one may 
obtain the most precise values by evaluating [A | and 
then inverting it. 
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Thermal Stress in Curved Beams (Continued from page 630) 


Here one may take 


V2? = (1/r) (d/dr) [r (d/dr)| (11a) 


Integration of Eq. (11) gives: 
= A, In (7 a) + Ad(r/a)? In (r/a) + A3(r a)? — 


Ea f nf Pr dr | dr (12) 


provided that a constant term which produces no stress 


isomitted. The boundary conditions 


= —(O Or) [(1/r) (06/00)] = atr=a,b (13) 


are automatically satisfied by this expression for ¢. 
The constants A;, As, and A; are to be found from the 
conditions that 


atr= a,b 


= (1/r) (dé/dr) = 0 


b b 
o,r dr = r(d*o/dr*) dr = 0 


It should be noted that the equilibrium condition 


(14) 


f ao, dr = O is automatically satisfied if o, = 0 
a 

atr = a,b. The expressions of Eqs. (5) for og and for 
the constants in question are readily obtained from the 
above conditions. An alternate derivation of Eqs. 
(5) may be obtained by superposition of the solution 
for a heated hollow cylinder and for a beam in pure 
bending, both of which may be found in reference 1; 
this fact was used as a check of the present work. 
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SUMMARY 


A study was made of the shock tube boundary layer between 
the initial shock wave and the contact region for shock pressure 
ratios Ps, of 2.75 and 8.0 with a corresponding Reynolds Number 
range of 0-50 X 10% and 0-100 X 10°, respectively. The Reyn- 
olds Number is based on the distance behind the shock wave 
and fluid properties evaluated at the wall. Shadowgraph, 
schlieren, and interferogram photographs and thin-film thermom- 
eter measurements indicated that the region was almost com- 
pletely occupied by turbulent boundary layer. Experimental 
density profiles and boundary-layer thickness measurements are 
presented together with calculated temperature and_ velocity 
profiles, displacement and momentum thicknesses, and local skin 
friction. Results are compared with the theory of Mirels in 
NACA TN 3712. In view of the assumptions of the theory the 
agreement is quite good. The velocity profiles at higher Reyn- 
olds Number for P2, = 2.75 follow a 1/5 power law closer than 
the 1/7 power law of TN 3712. For P2, = 8 all the experimental 
profiles are a 1/5 power. As a means of comparison the theory 
was redeveloped on the basis of the 1/5 power law for velocity. 


SYMBOLS 


Fig. 1 serves to define the flow regions in the shock tube. The 
notation Air/Air and He/Air denotes a gas combination across 
the diaphragm, the left-hand symbol representing the gas in the 
chamber and the right-hand symbol the gas in the channel. 


a@ = velocity of sound 
Cy» = local skin-friction coefficient 
C = contact surface moving to the right 


C, = specific heat at constant pressure 

C,, = wall material specific heat 

k, = wall material conductivity 

K = Gladstone-Dale constant = (index of refraction — 1)/p 

L = test-section span 

Lr = reduced test-section span due to boundary layer on test- 
section windows = L(1 — 26,) 


N= integer in exponent 1/N for the velocity profile of the 
turbulent boundary layer 

pressure 

R= gasconstant 

R- = backward facing rarefaction wave (particles enter from 


left) 


R, = Reynolds Number 

ry = recovery factor defined by Eq. (8) 

S = nondimensional fringe shift 

S = forward facing shock wave (particles enter from right) 
T = temperature 

t = time 
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u, v = velocities parallel to x, y axes 

u,0 = velocities parallel to j axes 

Uw» = velocity of wall in x, y coordinates 

w = velocity of initial shock wave in *, § coordinates 

X = distance from measuring station to diaphragm station 

x, = distance from initial shock wave to measuring station 

x, Y = coordinates with respect to shock wave 

¥, ¥ = coordinates with respect to wall 

z = coordinate normal to window with origin at the ray 
exit window 

6 = velocity-boundary-layer thickness 

6, = velocity-boundary-layer thickness on floor of shock 
tube 

6. = velocity-boundary-layer thickness on test section 
windows 

6* = boundary-layer displacement thickness defined by 
Eq. (10) 

oe = Prandtl Number 

@ = boundary-layer momentum thickness defined by Eq. 
(11) 

X = wavelength of light in angstrom units 

= kinematic viscosity 

p = density 

ps = wall material density 

7 = local shear stress 

Subscripts 
1 = undisturbed region ahead of the initial shock wave 


2, e = flow region behind the initial shock wave outside of the 
boundary layer 

m = quantity evaluated at a mean temperature 

w = quantity evaluated at the wall surface 


Dimensionless Quantities 


Ps, = pressure ratio across the initial shock wave, p»/p; 
Tx, = temperature ratio across initial shock wave, 72/7; 
MEASURING 
~ TION 
STATIO 
B=R 
= Te 
‘ 
(4) (3y(2 (1) 
R 
x 
vA 
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Fic. 1. The boundary layer in the shock tube. Note: 
subscript e is taken from reference 6 and subscripts 1 and 2 are 
taken from reference 8. 
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(b) WITH RESPECT TO WAVE 
Fic. 2. The coordinate systems used to study the shock tube 
boundary layer. 


INTRODUCTION 


I RECENT YEARS a great deal of interest has been 
shown in the shock tube boundary-layer problem. 
Much of this interest stems from the study of shock 
attenuation in the shock tube. The laminar-boundary- 
layer case was investigated theoretically in reference 1 in 
conjunction with a shock-attenuation study. In refer- 
ence 2 an attenuation theory was developed and com- 
parison with experiment indicated that the boundary 
layer was most likely to be of the turbulent type. 
Experiments carried out in reference 3 also agreed very 
closely with the theory of reference 2, however, both of 
these studies were concerned primarily with the atten- 
uation problem. The laminar boundary layer was 
investigated theoretically in reference 4, while refer- 
ence 5 presented a similar theory plus experimental 
verification. In reference 6 a_ turbulent-boundary- 
layer solution was added to the theory of reference 4. 
The present work is an experimental study of the 
turbulent boundary layer created on the floor of a 
shock tube by the passage of the initial shock wave 
down the tube; the results are compared with the 
theory of reference 6. This study was carried out in 
the UTIA 2 X 7 in. shock tube using a 9 in. field Mach- 
Zehnder interferometer to obtain boundary-layer 
density profiles. Profiles were measured in the turbu- 
lent boundary layer developed on the shock tube floor 
at various times after the shock passed a station located 
on the test section windows. The test region was 
bounded by the shock wave and the contact surface— 
ie., region 2 in Fig. 1. The boundary layer was 
studied behind shock waves of strength (P2,) 2.75 and 
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8.0. From the density profiles the temperature and 
velocity profiles were calculated, and, using these pro- 
files, it was possible to obtain displacement and momen- 


Boundary-layer thickness was meas- 
ured from the interferograms. Wall temperature in- 
crease was measured using a thin-film resistance 
thermometer. A more detailed description of these 
experiments is given in reference 7 and the pertinent 
shock tube theory will be found in reference 8. 


tum thicknesses. 


THE BOUNDARY LAYER IN THE SHOCK TUBE 


The shock tube boundary layer is somewhat different 
from the steady flow boundary layer found in the wind 
tunnel. However, in order to study this nonstationary 
boundary layer, a coordinate system has been adopted 
which essentially reduces the problem to one of steady 
flow, except that the wall has a finite velocity equal 
and opposite to the initial shock wave velocity. This 
method was adopted in references 1, 4, 5, and 6 and, 
since the present experiments were compared with the 
theory of reference 6, some of the notation of that re- 
port will be used to describe this boundary layer. 

Fig. 1 serves to define the regions in the shock tube 
and illustrate the shock tube boundary layer. Fig. 2 
presents the two coordinate systems used, one fixed 
with respect to the tube wall and the other attached 
to the shock wave. The two systems can be related as 
follows: 


FT=9 


x x — 


u 


(1) 


—w in the x, 


The wall moves with a velocity u, 
coordinate system. 

The velocity profiles in the x, y coordinate system 
corresponds to the usual profile found in two-dimen- 
sional steady flow. However, in this case, the choice 
of the coordinate system enables the boundary layer 
in the shock tube to be studied using proved steady- 
flow methods. Thus, in Fig. 2(b), the velocity profile 
obtained is complimentary to the profile of Fig. 2(a). 

The density and temperature profiles of the nonsta- 
tionary boundary layer are also opposite to those 
usually found for the wind tunnel case. In region 2 
outside the boundary layer the physical properties 
are po, py and T» (see Fig. 1). It was assumed in refer- 
ence 6 that the wall temperature remained constant 
at the room temperature, 7) (see also references 5 and 
9), which is also the temperature of the gas at rest 
ahead of the shock. The temperature increases from 
7T,, = T; at the wall (y = 0) to a value of 7, > 1 at the 
edge of the boundary layer (y = 6). Therefore, in the 
shock tube, the free-stream temperature is greater than 
the wall temperature, while in the wind tunnel the 
situation is usually reversed and the wall temperature 
exceeds the free-stream value. 

For the shock strengths of interest here, the gas may 
be assumed perfect with an equation of state p = pRT 
where R = constant, due to constant specific heats 
C, and C,. The assumption is also made that the 
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pressure in region 2 (pz) is constant across the boundary 
layer. The density then varies inversely as the tem- 


perature, and at the wall 
Pw/ P2 = = (2) 


Referring to the coordinate system of Fig. 2(b), 
the density ratio p, p2, is directly related to the ratio 
of the wall velocity and the velocity outside the bound- 
ary layer, “4, u,. This velocity ratio is used as a basic 
parameter in the theory of reference 6. The shock 
strengths of the present experiments correspond to 
values of wu, u, equal to 2.00 (Px = 2.75) and 3.50 
(Px = 8.00). 


EXPERIMENTAL METHOD 


Outline of Experiments 


The shock tube used for these experiments is de- 
scribed in reference 10. However, when the 9-in. field 
interferometer was built,'! a new test section’? was de- 
signed so that the interferometer could be used with 
this shock tube. It was the floor of this new test section 
that was used for the study of the turbulent boundary 
layer. 

Three series of experiments were carried out with 
Py = 2.75. At X = 77 in. (see Fig. 1), a He/Air 
combination was used which lowered the diaphragm 
pressure ratio but also reduced the test time. By 
using an Air/Air combination and increasing Y to 
141 in. the test time was tripled. However, the dia- 
phragm pressure ratio required was very close to the 
limits of the equipment. In both cases a channel pres- 
sure of 300 mm. Hg gave adequate fringe shift. The 
third series of experiments was at Y = 141 in. using 
a He Air combination. This meant a lower diaphragm 
pressure ratio, but in turn a shorter test time. Each 
of the above experiments was repeated as a check on 
repeatability. 

For a P2; of 8.00 it was necessary to use a He /Air 
combination to obtain this shock strength. The chan- 
nel pressure was 120 mm. Hg and the test section was 
located 141 in. from the diaphragm; the experiments 
were repeated twice. The pressure ratio of 8 was about 
the maximum attainable using a helium driver gas and 
sufficiently high channel density to provide adequate 
fringe shift with the 2-in. tube width. At both shock 
strengths the densities were high enough to give early 
transition from a laminar to a turbulent boundary 
layer. 

Records were obtained over the range of test time 
available by taking a series of spark photographs as 
the shock wave progressed downstream from the meas- 
uring station, producing a series of interferograms of 
the boundary layer in region 2. The diaphragm pres- 
sure ratio was adjusted so that the shock wave arrived 
at the measuring station with the required shock pres- 
sure ratio. One interferogram was obtained for each 
shock tube firing, the series being obtained by increas- 
ing the time delay controlling the spark. The runs 
were duplicated as closely as possible, so that the 
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SCIENCES 
records were obtained under near identical flow cond.) 
tions, as determined by the wave speed measured with | 
the light screens. Variations in wave speed were les 
than one per cent. 


Boundary-Layer Interferometry 


Various authors'*~'® have discussed the techniques 
of boundary-layer interferometry. 
effects may introduce errors into interferometric results 
The important effects for boundary-layer studies include 
misalignment of the test beam and the model, refraction 
of the test beam by the steep boundary-layer gradients 
and the presence of a boundary layer on the test section 
windows. The last is often referred to as a corner 
effect where the floor and window boundary layers | 
meet. 

A similar method to that used in reference 13 was | 
adopted here for aligning the test beam with the floor, © 
In the present experiments a special jig was con- 
structed with two plates at right angles to each other 
and perpendicular and parallel to the test section floor, 
In this manner, alignment was possible in both the 
vertical and horizontal planes by viewing the diffraction 
pattern about each plate and adjusting the interfer. | 
ometer until the pattern was symmetric. 
also served as a focusing device and a means of meas- 
uring the magnification factor in the focal plane (z = 
L/3). 

Because of the refraction and diffraction effects, the 
image of the floor of the test section on the photographic 
plate did not represent the true floor location with 
respect to the free-stream fringes. In order to locate 
the floor accurately, an inked reference line was placed 
on the test section windows, nearest the camera and 
several inches above the floor, well above the floor 
boundary layer. Before each set of experiments were 


In practice several | 


This jig | 


performed a no-flow calibration photograph was taken 
with the jig inserted in the test section. From this 
photograph a measurement of the distance from the | 
ink line to the floor was made, and the magnification 
factor determined from a measurement of the image 
of the jig. Having established the ink-line-to-floor 
measurement on the negative, this distance was meas- 
ured off from the line on each flow photograph of the 
set in order to locate the floor. Further details of this 
technique will be found in reference 7. 

Evaluation of some preliminary results indicated 
that any correction due to the refraction effects had a 
negligible influence in determining a density profile. 
However, it was also found that a correction was needed 
for the boundary layer on the test section windows. 
The density profile was evaluated by means of a rela- 
tionship developed in reference 16 which included 
the corner effects. This relationship for the present 
experiments is given below. 


= 1 + — Le) (1 — Tx) X 
+ N(y/s,) ] foy/(N + 1)] — 


The wall and window boundary-layer thicknesses were 
assumed to be identical (6. = 6,). 
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TURBULENT BOUNDARY LAVER 


Both horizontal and vertical fringe interferograms 
were taken during the preliminary studies. After ex- 
amining these two methods the technique of using 
closely spaced horizontal fringes, parallel to the test 
section floor, was adopted. One of the advantages of 
this method is that measurement in the y direction only 
is required for evaluation. These measurements were 
done on a comparator equipped with 10 power optics 
which could be read to 0.00005 in. The density pro- 
files were computed, using Eq. (3), on the University 
of Toronto digital computer FERUT. 


Wall Temperature, Transition, and Density Gradient 


In order to establish whether the theory of reference 
(, was applicable to these experiments, an experimental 
check was made on the assumptions of negligible stream- 
wise pressure gradient and wall temperature increase. 
The pressure gradient was estimated from interfer- 
ometer studies, while the wall temperature rise was 
estimated from thin-film resistance thermometer 
records. 

It can be shown" that such a thin film with its 
negligibly small heat capacity measures the surface 
temperature of the backing with negligible time lag. 
The thermometer was a gold film approximately a 
tenth micron thick, 1 8 in. wide by 3/4 in. long evapo- 
rated onto a backing of araldite plastic.* With the 
present films, the response time of the film to a shock 
wave was essentially the transit time of the shock wave 
across the 1/S-in. width. From the measured change 
in backing surface temperature, the temperature change 
of the steel wall could be computed. 

The thin-film thermometer was located 12 in. down- 
stream of the measuring station at the test section. 
The resultant signal was displayed on an oscilloscope 
and photographed; measurements were made from the 
record to determine the wall temperature rise. The 
transition from a laminar to a turbulent boundary layer 
was also indicated on these records. Shadow and 
schlieren photographs were taken as a means of check- 
ing on the transition point. 

The technique for determining free-stream density 
gradient entailed measurement of the fringe shift due 
to the density increase across the shock at a series of 
stations between the shock and contact front. Know- 
ing the density gradient, the corresponding pressure 
gradient could be estimated by assuming an adiabatic 
relationship between density and pressure. 

The fringe shift across the shock was measured in 
the first photograph of each set. However, some other 
reference was needed for the succeeding photographs. 
The reference used was the ink line drawn on the test 
section window. Before firing the shock tube, and 
with the channel held at the operating pressure, the 
zero-order fringe was aligned along the ink line using 
an AH3 steady source. Once this was done the cham- 
ber was pressurized, the plate holders inserted in the 
camera, the spark source replaced, and the tube fired. 

* The evaporation of the films was done by the Physics Depart- 
ment of the University of Toronto. 
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Schlieren photograph of the boundary-layer de- 
velopment behind the initial shock wave. The shock wave is at 
the extreme left (Ps, = 8.00, p; = 120 mm. Hg, He Air). 


PLATE 1. 


The time delay was adjusted so that a photograph 
was taken at a later time corresponding to another 
value of x, and the above procedure repeated. In this 
way the ink line represented the zero order fringe 
position before firing the shock tube. This reference 
could be used to measure the fringe shift for a number 
of stations behind the shock wave, and from this shift 
the density and pressure gradient could be determined. 


EXPERIMENTAL RESULTS 


Transition and Wall Temperature 


(a) Transition—Measurements taken from the schlie- 
ren and shadowgraph records indicated that only a 
very small portion of the floor boundary layer behind 
the shock wave was laminar. For Ps; = 2.75 with 
pi = 300 mm. Hg, transition started about 1.3 in. 
behind the shock but the boundary layer did not appear 
to become fully turbulent until after 4 to 6 in. With 
P» = Sand p; = 120 mm. Hg, transition started about 
0.75 in. behind the shock and at 3 in. a fully turbulent 
boundary layer had developed. Plate 1 is a schlieren 
photograph illustrating the growth of the boundary 
layer behind the shock. 

The thin-film thermometer records obtained on the 
7-in. side wall could also be used to estimate a transition 
point. Immediately behind the shock the heat transfer 
due to the laminar boundary layer varies as 1/+/1, 
which provides a step increase in wall temperature.” 
Upon transition to a turbulent boundary layer, the 
change in heat-transfer rate produces a marked change 
in slope on the thin film thermometer records. The 
lengths of side wall laminar boundary layers estimated 
in this way were 3 in. for P2, = 2.75, and lin. at P» = 8. 

The transition may be characterized by the Reynolds 
Number based on wall values 


Rery = (UeXr/ Ve) [1 — (Uy / Ue) |? (4) 


where xr is the distance behind the shock at which 
transition occurs, or by the same Reynolds Number 
based on free-stream values. 


Rer, = (UeXr/ te) |? (5) 

TABLE 

Method Po pi(mm. Hg) Rez, Rey, 
Optical 2.75 300 0.33 106 0.18 108 
(Tube Floor 

2” wide) 8 120 0.75 108 0.14 108 
Thin-Film 2.75 300 0.75 X 106 0.39 X 106 

Thermometer 
(7” side wall) 8 120 1x 10° 0.19 x 10° 
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The present transition results, in terms of Rer,, and Rer,, 
are summarized in Table 1. 

A systematic study of transition was considered 
outside the scope of the present work; however, the 
data obtained indicated that transition was completed 
about 4 in. behind the shock wave for the two shock 
strengths used. This was the location chosen for the 
first measuring station in both cases. 

(b) Wall Temperature Rise—In addition to informa- 
tion on transition, the film thermometer records per- 
mitted an estimate of maximum wall temperature rise, 
AT. It can be shown’ that the wall surface temper- 
ature rise is proportional to 4/1/8, where B, = 
Thus, the temperature rise A7\,, experienced by the 
steel shock tube walls in the present case, may be ex- 
pressed as 


AT, = ATr VBi/Bs (6) 


where subscripts a and s refer to araldite and steel, 
respectively. 

Measurements made at this laboratory gave a value 
of = 1.96 B.t.u. ft.-? (°F.-!) (hour)-"*. The 
corresponding value of 1/8, = 40.9 in the same units. 
Therefore, Eq. (6) becomes 


AT,,/ATp = 0.05 (7) 


The thin-film thermometer experiments may be sum- 
marized as follows: 


ATr aT, 
2.75 2.08°C. 0.10°C. 
8.00 0.56°C. 


This table shows that the wall remained at essentially 
constant temperature during these experiments. 


Streamwise Density Gradient 


Any ordered change in S from streamwise density 
gradient was masked by scatter. The maximum per- 
centage variation in pressure (assuming adiabatic 
pressure-density relation) corresponding to the 
maximum scatter in measured fringe shift is given in 
the summary below. 


Maximum Maximum 
Px Gas Combination % Var. in p % Var. in p 
2.75 He/Air 5 
8.00 He/Air 7 10 


This table shows that the scatter could mask a maxi- 
mum increase in pressure of 7 per cent for P2, = 2.75 
and 10 per cent for Px, = 8. In view of the results of 
various attenuation studies and experiments, some 
pressure gradient was to be expected. Pressure in- 
creases With increasing distance behind the shock have 
been reported by various authors—e.g., reference 2. 
The theory of reference 6 is for a shock advancing over 
a wall of infinite extent, whereas the flow in the present 
case is confined by the walls of the tube which causes 
a slight pressure gradient to be created. It is felt, 
however, that the possible pressure gradient existing 
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Density, Temperature, and Velocity Profiles 


Assuming constant pressure through the boundary 
layer, the ratio T/T, is equal to the density ratio 
p./p. Therefore, the measured density profiles are 
presented as temperature profiles in Figs. 3, 4, and 5, 
The conditions of each experiment are given on the 
respective Figures, each Figure representing the results 
of at least two sets of experiments. Due to the large 
temperature gradients existing at the wall it was not | 
possible to measure the wall value interferometrically, 
However, as shown by the thin-film thermometer re. 
sults, the wall temperature remains essentially constant 
at room temperature, which enabled the ratio 77, T, 
to be fixed. 

Figs. 6 and 7 are velocity profiles computed from the 
temperature profiles of Figs. 3 and 5, respectively, 
These profiles were calculated from the temperature 
results using the Crocco equation'’ in the following 
form 


T/T, = (T./T.) (1 + — | 
= (T,,/T,) [1 + — (8) 
c(6* 6)2/" (y 
where 
b = — 1, ¢ = [((T,/T.) — 


r(0) = Wom 


Although the Crocco relation was developed for a 
Prandtl Number of unity, the coefficients 6 and ¢ are 
presented in terms of a recovery temperature in keep- 
ing with the theory of reference 6, which assumes the 
Crocco equation should hold for values of the Prandtl 
Number not too far from one. 

Calculating the experimental velocity profiles for 
Prandtl Numbers other than one produced very little 
change in the results. Similarly, the theoretical tem- 
perature curves were virtually unchanged when calcu- 
lated for Prandtl Numbers of 0.69 (Px, = 2.75) and 
0.665 (Px = 8). These values were taken from refer- 
ence 24. 

The velocity profile may be expressed as a function 
of the similarity parameter (7 = y 6 or fr = (y 6*) 
(6* giving 


— tw) / (tte — ty) | = = 
for O = Cr (9) 


This relationship will be equal to unity for ¢7 > 1. 

In view of the lack of any existing experimental evi- 
dence for the turbulent boundary layer behind a moving 
shock wave, Mirels adopted the 1/7 power velocity 
profile from steady flow. This velocity relationship is 
plotted as a broken line on Figs. 6 and 7. The corre- 
sponding theoretical temperature profile is also plotted 
as a broken line on Figs. 3, 4, and 5. 

The experimental results of Figs. 3, 4, and 5 all tend 
to fall along a single curve over the range of measure- 
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The deviation from the 1/7 power is most noticeable 
on Figs. 4(b) (Px = 2.75) and 5 (Px = 8). On both 
of these curves there is very little scatter and the data 
tends toward a higher power profile. Noting this, an 
attempt was made to determine the experimental pro- 
file. A 1/5 power curve was found to fit the measured 
points more closely and is plotted as a solid line on Figs. 
3 through 7. Figs. 3 and 4(a) are similar in pattern, 
the experimental points falling between the theoretical 
curves, although the results of Fig. 4(a) do tend slightly 
toward the 1/5 power curve. In Figs. 4(b) and 5, 
however, there is very good agreement with the 1/5 
power profile with less scatter than the previous re- 
sults. Calculated velocity profiles are presented for 
each shock strength, and their behavior is nearly identi- 
cal to the measured temperature profiles. 

There may be some doubt as to the validity of the 
experimental velocity profiles. These profiles are 
deduced from density measurements on the assumption 
of constant pressure against the boundary layer using 
the Crocco equation [Eq. (8)]. This relation has not 
been proved valid for nonstationary flow. However, 
this relationship has seen widespread use for both 
laminar and turbulent boundary layers in steady flow. 
In reference 20 it was used in the same manner as used 
here—namely, to relate density, as measured inter- 
ferometrically, to velocity for a turbulent boundary 
layer. In reference 21 experimental verification of the 
Crocco relation was carried out by means of flight 
measurements of the boundary layer on a body of 
revolution. Therefore, there seems to be some justi- 
fication for its use here. 


Theoretical Consideration for N = 5 and 7 


As the experimental velocity and temperature pro- 
files seemed to follow a 1/5 power law, it was decided 
to carry out the method of reference 6 for VN = 5. The 
displacement and momentum thicknesses may be 
defined as follows: 


= [1 — (pu /p,u,) |dy (10) 


(pu p.u.) [1 — (u/u,) |dy (11) 
0 


By substituting Eqs. (8) and (9) into Eqs. (10) and (11), 
the following boundary-layer shape parameters are ob- 
tained: 


6 N(T, Tw) x 


(tw + [1 — te) Ty} (12) 
0/8 = N(T./Tz) (1 — (ttw/tte)] X 
(te ul + [1 = Up) 
[1 (Uw Ue) (13) 


where /y_,, /y, and /y+, are functions of b and c [see 
Eq. (8)] and are defined by the integral 
1 
i, = [Z*dz/ (1 + bZ — cZ?)| (14) 
0 


where a = N — 1, N, and N + 1. 
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Eqs. (8)—(14) were taken directly from reference (} 
except that the general power 1/N was used instead oj} 
the value of 1/7. In the remainder of this sectioy 
only the equations applying to V = 5 will be presented 
The corresponding equations for NV = 7 will be foung 
in reference 6. 

The integral form of the momentum equation is 


T d pu u d6 
Pelle? Pelle dx 
Assuming that the wall temperature is constant, 6 j 


can be considered independent of x and Eq. (15) may 
be written 


Tw/ Pete? = (0/5) (db/dx) (16 


In order to integrate Eq. (16), a relationship between 
7» and 6is required. The theory developed in reference 
6 uses the Blasius relationship for compressible turbu- 
lent flow over a semi-infinite flat plate, evaluating the 
fluid properties at a mean temperature, 7;,, as sug. 
gested by Eckert.”” 
The 1/7 power law, as given by Blasius, for the ve. | 
locity distribution in a smooth pipe is”* 


= 


In an attempt to follow the approach of reference 6 a 
curve was fitted to the famous results of Nikuradse for | 
smooth pipes, given in Fig. 20.4 of reference 23. There- 
fore, as the power, 1/5, was already established, it was 
necessary to determine the constant, C. A value for 
C of 6.2 gave a reasonable fit to the Nikuradse results | 
similar to that for C = 8.74. This value resulted in | 
the following expression for the skin friction: 


(17) | 


pu? = 0.0488 (18) 
For a shock tube, this equation becomes 


Tw/ Pele? = 0.0488 — X 
— (Uy/Ue) (v./u6)'/* (19) 


Substituting Eq. (19) into Eq. (16) gives 


= 0.130 {¢[1 — (wu Xx 
[1 — (tty /tte) (20) 


Combining Eqs. (19) and (20) with Eq. (4a), where 47 
is replaced by the more general value x, results in the 
local skin friction 


R, = 0.195 
i, 6 t — 1% 
Theoretical curves of 6, 6*, 6, and c,,, are shown on Figs. 
8, 9, and 10 for both the 1 5 and 1/7 power law profiles. 


where 


Boundary-Layer Thickness, Displacement, and 
Momentum Thickness 


Figs. 8 and 9 present plots of 5, 6*, and @ versus x, and 
Re,, {based on wall values, Eq. (4a)]. The boundary- 
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Fic. 8. Theoretical and experimental variation of the boundary 
thickness 6, displacement thickness 6*, and momentum thickness 
4, with the distance behind the shock wave x, and the Reynolds 
Number R.,,. Pa = 2.75, pr = 300 mm. Hg, tw/u, = 2.00. 
Note: the Reynolds Number in Figs. 8, 9, and 10 are evaluated 
for the fluid properties at the shock tube wall. 


layer thickness, 6, came out of the program for deter- 
mining density profiles, this thickness being arbitrarily 
defined as the value of 6 corresponding to a density ratio 
of 0.996. Experimental values of 6* and @ were deter- 
mined from the density and velocity profiles by numer- 


' ical integration of Eqs. (10) and (11) using FERUT. 


The theoretical and experimental values of 6* and @ 
agree reasonably well for Ps; = 2.75. By the definition 
of 6* and @ these values are both negative and in these 
coordinates @ is larger in magnitude than 6*. The 
theoretical curve for NV = 7 matches the experimental 
data somewhat better than for V = 5. 

For Ps, = 8.00, the two theoretical curves differ 
but slightly and, particularly for 0, both fall below the 
experimental results. The reason for this is not known. 
The large scatter around x, = 28 in. for 6* and @ 
is probably due to the conditions that exist in the 
shock tube at this value of x,. This station is 
very close to the contact surface or cold front. The 
region behind the contact surface is not uniform and 
it is possible that disturbances are transmitted forward 
through the boundary layer to perturb the flow in re- 
gion 2. 

Too much significance cannot be attached to the 
comparison of theory with experiment for boundary- 
layer thickness in Figs. 8 and 9, as the experimental 
thickness has been arbitrarily taken as the value of y 
for pp p = 0.996. The theory for V = 7 overestimates, 
and that for N = 5 somewhat underestimates, this 
experimental thickness. 


Skin Friction 


Fig. 10 presents the results for local skin friction 
plotted against Reynolds Number. This Reynolds 
Number is identical to that defined by Eq. (4a). The 
experimental points were obtained by measuring the 
slope of a mean curve through the experimental @ versus 
\, results in Figs. S and 9 and then using the following 
equation for cy: 


Cyr,, = (d0/dx) (1/2) pu(tte — Ue)? | (22) 


HIND THE INITIAL SHOCK WAVE § 651 


The values of cy, are negative in the present coordinate 
system. 

The agreement between theory and experiment is 
fairly good for Px, = 2.75 [Fig. 10(a) |, the experimental 
skin-friction values falling between the two theoretical 
curves. However, for P:, = 8 the experimental points 
lie well above the two theoretical curves. This is to be 
expected in view of the corresponding results for mo- 
mentum thickness shown in Fig. 10(b), which indicates 
a greater slope than the theoretical curves. The last 
group of points at x, = 28 in. were not used in deter- 
mining the slope for the skin-friction calculation. 
The Reynolds Number was not corrected for the small 
laminar portion of the boundary layer in Figs. 8, 9, 
and 10. However, the correction will be small in rela- 
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tion to the very large Reynolds Number obtained as 
defined by Eq. (4a). 


CONCLUSIONS 


The assumptions of reference 6 were as follows: 
(1) zero pressure gradient, (2) velocity profile to follow 
a 1 7 power law, (3) Croeco relation related density 
or temperature to velocity, (4) wall temperature re- 
mained constant at room temperature, and (5) Blasius 
relation for the shear stress was extended to the moving 
wall case. 

In the present experiments, a small streamwise pres- 
sure gradient probably existed with pressure increasing 
with distance behind the shock wave. An attempt 
was made to measure the streamwise gradient by inter- 
ferometry, but any ordered change was masked by 
scatter. Measurement of the steel wall temperature 
rise by means of a thin-film resistance thermometer 
showed it to be negligibly small. 

In general, the theory of reference 6 compares favor- 
ably with the experimental results obtained, consider- 
ing the assumptions of the theory and those inherent 
in deducing velocity from the experimental density 
profiles. In particular, good agreement is obtained 
for the weaker shock case with ?:; = 2.75. For the 
stronger shock case with P2, = 58.00, experimental 
displacement and momentum thicknesses and skin 
friction are found to be greater in absolute value than 
predicted by the theory. 

At the higher Reynolds Numbers, and particularly 
for Px = 8.00, the experimental temperature and 
velocity profiles followed a 1/5 power law more closely 
than a 1 7 powerlaw. The method of reference 6 was, 
therefore, applied for a 1/5 power law to calculate 
displacement and momentum thicknesses and_ skin 
friction. The calculated results for the 1/5 power law 
differ very little from those for the 1 7 power velocity 
profile and show little, if any, improvement in agreement 
with the experimental results. 
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Upper Atmosphere Atomic-Oxygen Power 
Plant 


Sterge T. Demetriades 


Technical Consultant, Aerojet-General Corporation, 
Azusa, Calif. 


May 15, 1958 


| Shgerpremmense INTEREST has been generated about the family 
of ‘‘excraft”’ 

the Astronautics Laboratory of Corporation 
under sponsorship of the Directorate of Advanced Studies of the 
Air Force Office of Scientific Research. Recombination of the 
atomic oxygen available in the upper atmosphere* has been con- 


(exospheric aircraft) that is being studied at 
Aerojet-General 


sidered as the source of energy for excraft propulsion. 

The author and coworkers ? have presented theorctical esti- 
mates of the recombination kinetics of atomic oxygen. These 
results can be summarized by a homogeneous (volume) recom- 


bination rate given by 


[d(Oz) = Ri(O)*(u) {1 /Kea(O)2|{ + 


[(O3)?/(O2) |, moles/sec.-ec. (1) 


where (Q), (Oz), (Os), and (#) are the concentrations of oxygen 
oxygen molecules, and ozone molecules, and the total 
; ky is the three-body re- 
combination rate constant; k» is the rate constant for ozone de- 
is the equilib- 
is the equilibrium constant 


atoms, 
concentration respectively in moles/cc 


composition by the mechanism O; + O — 20,; 
rium constant (Oz) /(O)?; and As 
(Oz) (O)/(O;). These constants are given by 


* 1 prefer the term exosphere for the lower upper atmosphere and the 
term phygosphere (Greek phygé = escape) for the extreme upper atmos 
phere where the mean free path is larger than the scale height 
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ec.2/mole?-sec., 300°K. < 
T < 4,000°K. (2) 


ky = 3 x 10 +1) (300/T)*, 


ke = 1.07 X 10°38 (7/500) 10-7), ee. /mole-see. | (3) 
Kea = (Oz) /(O)2 = 101(%,700/ 7) ~1. 390] 


and 
Ke = (O) /(Oz) = 7) +6.93 logs, RT} (4) 

where R = 82.06 atm.-cc./°K.-mole; and T = °K. The un- 
certainty of (+1) on the exponent of Eq. (2) arises mostly from 
the different efficiences of various possible third bodies. The 
effect of nitrogen and nitrogen oxides on the recombination rate 
of atomic oxygen is still uncertain. 

The Demetriades-Farber estimates of heterogeneous (catalytic) 
recombination on surfaces can be summarized by 


[d(O>») /dt|yr = 1.5 X 105(O)24/T + 
9 X 1027” O)fYT, 


The efficiency f ranges from 0.001 for glass or aluminum to about 
1.0 for the noble metals and some oxides. 


moles/sec.-em.? (5) 


TABLE 1 


Summary of Upper Atmosphere Atomic Oxygen Recombination 
Power-Plant Spectrum 


Velocity Subor bital Orbital 
Moch Number 
at 00 km _|0, 2, 4,6, 8, 10, 2, If, 16, iB, 20, 22, 24, 3, 28, 30, 32 
Catalytic 
Type of Catalytic with ram comp 
Power Plant Pure 
++——Thrust at IOO km less than 40 dynes/cm* of inlet-———— 
compression 
Lift too low + 
(L<0.01 gm/cm?) (L<0.08 gm/cr? ) 
»-—Compression required supersonic throughout poker plant— 
Difficult inlet ond diffuser design 
*—Skin friction too high with present estimates 
————Wove or pressure drag losses could be severe ———~ 
Length of pure compression power plant 
becomes excessive with present estimates 
losses from 
catalytic walls 


“Lift per cm? of |ift-producing surtace 


The recombination power-plant studies conducted at the Astro- 
nautics Laboratory and now in process of publication are sum- 
marized in Table 1. Considerable uncertainty exists about the 
effect of surface recombination on the skin friction, wave drag 
at 100 km., diffuser design (compression ratios of 1,000—5,000 are 
required), lift, availability of extremely light structures, etc. I 
doubt, however, if these parameters can be improved by the re- 
quired order of magnitude soon. 

In the atomic oxygen layer one might be able to extract lift 
from a foil at zero angle of attack coated with a catalyst (of low 
accommodation coefficient) on one side and a noncatalyst (of 
high accommodation coefficient) on the other, and thrust from 
a “wavy” foil coated with this catalyst on the forward-facing 
surfaces of its waves and a noncatalyst with high accommodation 
coefficient on the backward-facing surfaces. 
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Static Friction in Dynamic Systems 


Carl F. Branson 
Senior Aerodynamics Engineer, Lockheed Aircraft Corporation, 
Georgia Division, Marietta, Ga. 


May 19, 1958 


opm USUAL CONCEPT of static or coulomb friction is a constant 
force with its sign opposite to the velocity. The explanation 
of some self-excited vibrations as given in reference 1 requires 
the force to decrease with the velocity, retaining the opposite 
sign, however. Such definitions are sufficient for predicting the 
effect on single-degree-of-freedom systems, except for the time 
history. 

In determining the damping and response of a system with 
two or more degrees of freedom, the behavior of the friction when 
the velocity is zero becomes important. Consider the friction 
acting in one degree of freedom only of a two-degree-of-freedom 
system. The friction, when the velocity is zero, is exactly equal 
in magnitude and opposite in direction to the accelerating forces 
until the maximum friction is exceeded. Once the maximum 
friction is exceeded, an acceleration takes place and the velocity 
The usual definition of friction then 
The accelerating 


becomes other than zero. 
applies until the velocity becomes zero again. 
forces referred to are the spring, cross coupling, and exterior 
force, if any. 

The equations for the simple two-degree-of-freedom system 
can be written: 


+ F + + = 0 
Mois + + = 0 


where a; and a» are the coupling terms. The effect of the friction 
can be determined from the energy exchange as shown in phase- 
plane plots. The plot with no friction is shown in Fig. 1. The 
relationship between x; and xe is simply a straight line of slope 
dx;/dx. = a;/k;. If we consider a2 the primary variable, the 
addition of x; as another degree of freedom changes the spring 
rate but has no effect on the damping. With the addition of 
static friction to x,, however, the motion can either add energy 
to or subtract energy from the x. system, respectively, either de- 
creasing or increasing the damping, depending on the relative 
signs of a, and ae. 

The phase-plane plots of Fig. 2 are drawn for very smail inertia, 
m,—.e., 


m, = 0 


in order to more clearly show the energy relationships due to fric- 
tion. The energy added to the x2 system would be 


AE = f 
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where a resultant positive sign for AE would indicate energy 
added and decreased damping. Consider the case where a; and 
a are both negative. From A to B, energy is added to the sys- 
ten. From B to D energy is subtracted in a larger amount 
(since dx. is negative). The resultant effect of the friction is 
damping, the other half cycle giving identical effect. The oppo- 
site would be true if the sign of a2 were changed to positive. Now 
consider a, and ay both positive. From A to E energy is again 
added to the system and from E to D energy is subtracted with a 
resulting increase in damping. It is obvious that positive a 
and negative a2 would give the reverse. Any addition of energy 
must come from some outside source such as the air stream in the 
case of airplane dynamics. 

It is interesting to note the part the cross coupling terms 
play. When the signs are alike, the effect of friction is increased 
damping, and when they are not alike, negative damping occurs. 
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Vaporization Into a Hypersonic Laminar 
Boundary Layer 


Sinclaire M. Scala* 
Research Engineer, Missile and Ordnance Systems Department. 
General Electric Company, Philadelphia, Pa. 


May 27, 1958 


SYMBOLS 


a@ = vaporization coefficient 
cp = > nCp, = frozen gas specific heat 
n 
Cn = mass fraction of species ” 
YD, = diffusion coefficient of species 1 
jn = pnb diffusion flux vector of species 1 
Le = Lewis Number 
mi = interphase mass transfer 
M, = molecular weight of species » 
p = static pressure 
Pi = Prandt! Number 
R = universal gas constant 


* The author wishes to acknowledge the many helpful discussions held with 
Dr. Guido Vidale of the Aerosciences Laboratory, M.O.S.D., General Elec 
tric Company, Philadelphia, Pa ca 

This analysis is based upon work performed under the auspices of the 
USAF-B.M.D. contract No. AF 04(645)-24 
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T = temperature 
Un = absolute velocity of species » 
(due/dx)s = stagnation point velocity gradient 
6 = T/T, = dimensionless temperature 
p = density 
“ = viscosity 
” = similarity variable 
Subscripts 
A = atom 
e = edge of boundary layer 
i = interface between gaseous and condensed phases 
eq = equilibrium 
M = air molecule 
k = vaporizing species 
n = denotes differentiation with respect to n 


Fata of surface melting or sublimation requires a detailed 
consideration of the vaporization processes at the interface 
between the condensed phase and the gas phase. In Fig. 1 are 
indicated the various physicochemical processes which occur in 
the dissociated hypersonic boundary layer, in the mean free path, 
(shown to exaggerated scale) near a surface undergoing vaporiza- 
tion. 

Heterogeneous atom recombination has been considered pre- 
viously,! when mass transfer was absent. However, when 
vaporization occurs, a new set of constraints is applicable. Al- 
though these boundary conditions were utilized recently? in 
considering the problem of a melting surface, few details were 
given due to space considerations. Because of the current in- 
terest in this problem, the surface constraints will be presented 
here. In particular, the relationship between the interphase 
mass transfer m;, the diffusion rate of the vaporizing species 
je, the equilibrium vapor pressure px;,,, and the partial pressure 
px; will be discussed. 

A body immersed in a high-temperature, high-velocity stream 
behaves initially as a pure heat sink. When the surface reaches 
a critical temperature, selected particles at favorable sites will be- 
come detached from the surface lattice structure. During the 
microscopic collision processes which follow, a number of par 
ticles are readsorbed by the surface, while others are transported 
away by convection and diffusion in the fluid stream. 

If one assumes that the forward rate of vaporization from the 
surface is independent of the partial pressure in the gas phase 
pxi, and that the rate of condensation is proportional to the partial 
pressure px;, then kinetic theory yields: 


(a Me /V 24 (Pig, — pki) (1) 


as the net rate of vaporization where px; is evaluated within one 
mean free path from the surface. This equation states that the 
actual partial pressure px; is less than the equilibrium vapor pres- 
sure at the interface temperature 7;, when a finite interphase 
mass transfer occurs. 

The interphase mass transfer, by convection and diffusion in 
the gas phase may be written 


= = (pava + + (2) 
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NORMALIZED INTERPHASE MASS TRANSFER 
Fic. 2. Mass transfer of vaporizing species in the hypersonic 
laminar boundary layer. Stagnation point altitude 100,000 ft., 
Mach Number 20. 


Heterogeneous atom recombination requires 
(pava)i = (3) 
which yields the simple result: 
mi = (pits (4) 
Noting that Din = (), rearrangement of either Eq. (3) or Eq. (4) 
vields the constraint : 
mi; = — cri)! (5) 


Upon introducing Fick's law which is the first approximation to 
the diffusion flux of the kth species, 


Jk = (6) 
Eq. (5) may also be written in terms of the similarity variables 
(3) 


V — cr;)| (7) 


at the stagnation point of an axially symmetric body. In the 
steady state, Eqs. (1) and (5) are equal. Hence, if the equilib- 
rium vapor pressure can be determined as a function of tempera- 
ture, for example, by integrating the Clausius-Clapeyron equa- 
tion for the particular vaporizing material, and if the vaporiza- 
tion coefficient a of the surface is known, then and only then may 
one compute the actual partial pressure px; of the vaporizing 
species at the interface, a posteriori. If one merely assumes that 
the equilibrium vapor pressure is equal to the partial pressure, 
one may actually overestimate the vaporization rate, while 
solving the boundary-layer equations. The discussion presented 
above should be critically compared with the recent conclusions of 
Bauer and Zlotnick.* 

Because of the similarity between the energy and diffusion 
equations, for a frozen hypersonic boundary layer, it can be shown 
that :3 


(Cpi/Cpe) On; — 1) (8) 


Since the author has already correlated the functions 6; and 
(€pi/Epe) as a function of flight speed, altitude, mass-transfer rate, 
and surface temperature, for the binary mixture, hypersonic 
laminar boundary layer with mass transfer,> Eq. (8), when in- 
serted into eq. (7), yields the typical results shown in Fig. 2 
(which appeared previously as Fig. 4 of reference 2). The com- 
plete correlation of the function m; as a function of c,; for arbi- 
trary flight conditions is presented in reference 6 and may be 
utilized in the analysis of arbitrary materials, vaporizing under 
the particular environmental conditions of interest. 


REFERENCES 


'Seala, S. M., Hypersonic Heat Transfer to Catalytic Surfaces, Journal of 
the Aeronautical Sciences, Vol. 25, No. 4, pp. 273, 274, April, 1958. 

2? Scala, S. M., and Sutton, G. W., The Two Phase Hypersonic Laminar 
Boundary Layer—A Study of Surface Melting, presented at the 1958 Heat 


OCTOBER, 1958 


Transfer and Fluid Mechanics Institute Symposium, pp. 231-240, June 
1958. 

3 Scala, S. M., Vaporisation into a Hypersonic Laminar Boundary Layer 
General Electric Co., Aerophysics T. M. 84, January, 1958. 

‘ Bauer, Ernest, and Zlotnick, Martin, Evaporation into a Bounday, 
Layer, AVCO RAD Report, May, 1958 

5 Seala, S. M., Transpiration Cooling in the Hypersonic Laminar Bounj 
ary Layer (to be published). 

Seala, S. M., and Vidale, G., Vaporization Processes in the Hy person, 
Laminar Boundary Layer, General Electric Co., Aerophysics Research Mem 
orandum No. 14, July, 1958. 


On Simplified Fuselage-Structure Stress 
Distributions 


W. R. Jensen 

Structural Methods Engineer, Grumman Aircraft Engineering 
Corporation, Bethpage, Long Island, N.Y. 

April 3, 1958 


id IS GENERALLY AGREED that, in order to obtain an accurate 
stress distribution in a semimonocoque structure, a rather 
complex idealized structure is required. This is no longer an in- 
surmountable obstacle to stress analysts, thanks to the large 
digital computers now available. Many engineers doing this 
work nevertheless believe that such problems should still be kept 
as simple as possible, avoiding in this way excessive preparation 
of data, increased possibility of making errors, and voluminous 
answers. However, little information is available as to what con- 
stitutes the minimum in an adequate idealization. The present 
note attempts to shed some light on this matter. 

A well-known technique employed in idealizing these semi- 
monocoque structures consists of combining the stringers and 
skin into a smaller number of equivalent stringers and combining 
the frames into a smaller number of equivalent frames. The skin 
in between is assumed to be a membrane which can carry shear 
only. In an attempt to decide how many equivalent stringers 
should be used, several such idealizations of an NACA test cyl- 
inder! were made, and solutions were obtained analytically. A 
comparison of these results with a more comprehensive analysis 
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ind with test data! indicates that a surprisingly small number of 
equivalent stringers may be adequate for the design of most 
fuselages 

The actual structure analyzed (see Fig. 1) consisted of a 
radially loaded ring supported by a four-bay cantilever shell. The 
shell had no stringers and the idealization consisted of considering 
the longitudinal load carrying material to be concentrated at six 
or eight symmetrically placed locations at the skin line (see Fig. 
1). The more comprehensive analysis! * considers the longitu- 
dinal bending material to be continuously distributed around the 
shell, This latter analysis thus shows what can be expected 
when the number of stringers of the idealized structure is increased 
indefinitely. 

If the loaded ring of the structure were rigid, the stress dis- 
tribution in the shell would be according to engineering beam 
theory. Deformation of the loaded end ring induces self balanc- 
ing loads in the shell. These redundant distributions can be 
found by using a least-work method? and superimposed on the engi- 
neering beam-theory distribution. Bending energy of the loaded 
ring and reinforcing rings, shear energy of the skin, and energy 
due to the longitudinal bending stresses (stringer stresses) were 
considered. 

The six- and eight-stringer results are compared graphically 
with the more comprehensive solution, with test data, and with 
engineering beam-theory distributions in Figs. 2, 3, and 4. The 
ring bending moments for the six-stringer, eight-stringer, and dis- 
tributed-stringer analyses are so close that they can only be 
plotted as single ctftves (see Fig. 2). The test data are in excellent 
agreement. 

In Fig. 3 the curves represent the longitudinal bending stress 
for the distributed-stringer analysis. The stresses for the six- and 
eight-stringer analyses are plotted at their proper locations; as 
can be seen, faired curves through these points would agree very 
well with the distributed-stringer result. There are few test data 
for this comparison, however, excellent agreement is shown for 
the data available. 

The shear-stress distributions in the shell are shown in Fig. 4. 
The six- and eight-stringer cases are step curves as a consequence 
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Fic. 2. Ring bending moment. 
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Fic. 8. Shell longitudinal bending stress. 
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of the idealization. Within this limitation, the agreement is 
satisfactory. 

It can be concluded that, at least for a structure of the propor- 
tions and type of loading considered here, a very much simplified 
idealized structure with a small number of equivalent stringers 
will give satisfactory design information. 


REFERENCES 


! Kuhn, Paul, Duberg, John E., and Griffith, George E., 7he Effect of 
Concentrated Loads on Flexible Rings in Circular Shells, NACA, ARR No 
L5H23, 1945 

2 Wehle, L. B., and Lansing, Warner, A Method for Reducing the Analysis 
of Complex Redundant Structures to a Routine Procedure, Journal of the Aero 
nautical Sciences, Vol. 19, No. 10, pp. 677-684, October, 1952 

3 Hoff, N. J., Stresses in a Reinforced Monocoque Cylinder Under Concen 
trated Transverse Loads, J. Appl. Mech., Vol. 66, pp. A-235-A-239, Dec., 1944 


Further Remarks on the Torsional Rigidity of 
Thermally Stressed Wings 


R. L. Bisplinghof 
Professor of Aeronautical Engineering, Massachusetts Institute of 
Technology, Cambridge, Mass. 


June 5, 1958 


tae AUTHORS—e.g., references 1 through 6—have pointed 
up the important influence of finite deformations and initial 
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imperfections on the torsional rigidity of thermally stressed 
wings. The present brief note is intended primarily to com- 
municate a comparison of theory with experiment for finitely 
twisted flat plates subjected to a chordwise temperature distri- 
bution. Referring to Fig. 1, suppose that a flat plate of thickness 
h, chord 26, and length 2a (a >>b) with a chordwise distribution 
of temperature described by the even function 7(x2) undergoes 


finite twisting deformations due to a torque, J/;. It is known 
that when warping restraint is neglected, the torque is related 


to the twist per unit length, 6, by 
+b 
M, = — | (1) 


where D = Eh*/i2(1 — vy?) is the flexural rigidity, y is Poisson’s 
ratio, E is Young’s modulus and N;,;(x2) is the stress resultant in 
the x, direction. It can be shown that the stress resultant has 


the simple form? 
( Eh6?/6 (3x2? — b?) + 
b 
(a@Eh/2b) Lf, T(x2)dx. — (2) 


where a@ is the coefficient of thermal expansion. Assuming 

T(x2) = Ty + ATf(x2), we obtain by substituting Eq. (2) into 

Eq. (1), 

[A1,/4bD(1 — = 1 — (AT/AT) + (4/15)(1 + )(6402/h?) 
(3) 

where AJ is the temperature differential required to produce 

thermal buckling. 


b b 
ar = +1) [ — (04/3) | 


(4) 


The relative influence of the temperature differential, A7, and 
twist 1ate, 0, is indicated by the right-hand side of Eq. (3). 

In Fig. 2, there is shown a comparison of Eq. (3) with experi- 
ments? in which a 7075 T-6 aluminum alloy flat plate is twisted 
under the three conditions of AT/AT = 0, 0.4 and 1. A chord- 
wise temperature distribution was produced by using identical 
line heat sources along the edges and a line heat sink along the 
mid-chord. The temperature distribution obtained in this way 
is shown by the sketch of Fig. 1. 

In general, close agreement between Eq. (3) and the twisting 
experiments was found provided the initial imperfections were es- 
sentially eliminated and provided the thermal buckling tempera- 
ture differential, AT, was obtained with precision. The com- 
parison shown in Fig. 2 represents results obtained by giving care- 
ful attention to these two factors. Although reasonable agree- 
ment was found to exist between Eq. (4) and the experimentally 
obtained value of AT, the latter was used in preparing the solid 
lines in Fig. 2. 


Fic. 1. Experimental model. 
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Polynomial Expressions for the Specific Heat 
and Prandtl Number of Air 


Donald W. Boekemeier 
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May 22, 1958 


i pow PURPOSE of this note is to present polynomial expressions 
and their associated coefficients for the specific heat and the 
third root of the Prandtl Number of air as functions of tempera- 
ture. 


SYMBOLS 


Cp = specific heat of air 
D = coefficients in polynomial expression of specific heat 
E = coefficients in polynomial expression of (Pr)'/* 
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READERS’ FORUM 


Temp. Range = _300°-1200°R 


Pressure = 0.10 AIM 


+3.418923 
+l,.8331520 x 1 
-1.6621892 x 10-6 


-1.322860 x 10712 


Table I 


Temp. = 1200°- 


Pressure = 0.10, 1,00 and 4.00 ATM 


+3.221,0709 
+2.51,5360, x 1074 
+2,5222882 x 1077 
2013920 x 10-10 
+1.700360 x 10714 


on 


Coefficients for Polynomial Expression of Specific Heat 


Temp. Range = 3000°-51,00°R 
Pressure = 0.10 AIM 


= +6,9192916 x 
-9.3467711 1072 
= x 10-5 
= -1.1268396 x 10-8 
= +1,8)87752 x 10-12 
-8.746370) x 10717 


D3 = +264380548 x 10-9 


+2,1;858005 x 10-16 
Error = 


Pressure = 1.00 ATM 

Do * +36 71 

Dl = -2.0430912 x 10-4 
Do = x 10-7 

D3 = +8,144198) x 10-10 
Di, = -429535650 x 10-13 
Ds = +9.0255235 x 10-17 
Maximum Error = 0,08% 


io = 

Ih = -2.865688) x 10-3 

D2 = +5.8753994 x 10-6 

D3 = -529300099 x 10-9 

Dj, = +3.1877628 x 10-12 
D5 


Ds = 0,0000000 
Maximum Error = 0.50% 


Di = -3.8410018 x 10-2 
De = +2.133168) x 10-5 
D3 = -5.689,998 x 10-9 
Dj, = +7.2132357 x 10-13 
D5 = -3.3329913 x 10-17 
Maximum Error = 0.70% 


Pressure = 200 AT™ 

D1 = -2.54,33536 x 1072 
Do = +1.3943727 x 10-5 
D3 -3.662),308 x 10-9 
Di, = +.5962131 x 10-13 
Dg = -2.12998)5 x 10-17 

Error = 0.35% 


(Pr)'* = Prandtl Number to the one-third power 
R = universal gas constant, 1716.5 ft. Ib.¢/slug °R. 
= temperature, °R. 


Since many heat-transfer problems are presently being solved 
by electronic computing machines, it is necessary that the ther- 
modynamic and transport properties of air be either stored in the 
machine or calculated from equations. Properties such as the 
thermal conductivity and viscosity have accepted expressions as 
a function of temperature, for instance, reference 1. The spe- 
cific heat and Prandtl Number do not have expressions commonly 
given; they are usually just tabulated. 

Before the coefficients for these two expressions are determined, 
it is necessary to resolve what are the more widely accepted and 
agreed upon data presently available. After reviewing refer- 
ences 1 and 2, it is found that reference 1 accounts for the disso- 
ciation of the air at the higher temperatures whereas reference 2 
does not. All other investigation verifies the fact that dissoci- 
ation has a major effect on the properties of air at high temper- 
atures and that the data given in reference 1 are the best informa- 
tion presently available in the temperature range 200 to 5,400°R.; 
therefore, these data are used for determination of the coefficients. 

Three different constant pressures are taken for the determi- 
nation of the coefficients for the specific-heat polynomial expres- 
sion, the pressures being 0.10, 1.00 and 4.00 atm. The data are 
divided into three temperature ranges, and individual polynomial 
expressions are written for each pressure in these ranges. The 
ranges in which these data are divided are: 300—-1,200°R., 
1,200-3,000°R, and 3,000-5,400°R. All three pressures have 
nearly the same value of specific heat at a given temperature in 
the middle region, therefore, one set of coefficients suffices for the 
three. 

The polynomial expression taken is of the general form: 


C,/R = Do + DiT + DoT? + D3T? + + 


The constant coefficients for each pressure and temperature 
range are listed in Table 1 with the maximum deviation error 


Coefficients for Polynomial Expression of Prandtl Nunber 
Temp. Range = 300°-1800°R 


Pressure = 1,00 ATi 

Eo = +0.93790000 

Ey) -0.95038890 x 1074 
E2 -081,9383310 x 1077 
= +0.16898148 x 1079 
= -0.10339506 x 10-12 
BS = +0,20233196 x 10716 
Maxirmn Error = 0.10% 


indicated. It must be remarked here that even though the 
coefficients listed are for a certain temperature range, this does 
not sharply limit the use of these coefficients to this region. Their 
usage may be extended, but as it is extended the deviation error 
becomes larger, so it is advisable to limit this extension to 100°. 

Data for the Prandtl Number from reference 1 are given only 
for atmospheric pressure over the temperature range 300°—1,800° 
R. Only one set of coefficients is used for the entire range. 

The polynomial expression taken is of the general form: 

(Pr)'” = Ey + + + EsT* + E,T* + 

The constant coefficients for the given pressure and temperature 
range are listed in Table 2 with the maximum error indicated. 
It might be noted that the indicated maximum error is less than 
the accuracy of the tabulated data source. 
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On a Complementary Energy Principle in 
Linear Thermoelasticity* 


G. Herrmann 

Associate Professor of Civil Engineering, 
Structures, Columbia University, N.Y 

June 9, 1958 


Institute of Flight 


= EQUATIONS of linear thermoelasticity, using in part the 
formulation given recently by Biot,! consist essentially of 
the following set of four equations 


— = 0 (1) 
0; + pTs:/k = 0 (2) 
(ui; + = (3) 
Sic + = —Wo (4) 


Use is made of the summation convention. A comma preceding 
a subscript indicates differentiation with respect to that variable, 
except that 7 (or j) after a comma stands for differentiation with 
respect to the Cartesian coordinate x; (or x;). 

In the above equations, 7;; = 7); denotes the components of 
stress if the temperature increment @ above a reference absolute 
The total stress, in general, is thus 
Further, “,; and s; indi- 


temperature 7 vanishes. 
Tij — B06;; (6;; is the Kronecker delta). 
cate the components of solid and entropy displacement, respec- 
tively, & is the coefficient of heat conduction, and p = d/dt is a 
time operator. W is Biot’s thermoelastic potential, which con- 
sists of the sum of the isothermic mechanical energy (expressed 
here as a quadratic function of the stresses 7,;) and the term 
c02/2T, where c is the specific heat of the unit volume in the 
absence of elastic deformation. It is recalled that the entropy 
displacement s; is related to entropy s by the equation s = —s;,;. 

Eq. (1) represents the stress equation of equilibrium, in which 
the body force was omitted for simplicity. Eq. (2) is a particular 
form of the heat conduction equation. Eq. (3) is the (isothermal) 
stress-displacement relation (Hooke’s law), and Eq. (4) may be 
called the ‘‘thermal’”’ stress-displacement relation. 

Biot has shown in reference 1 that Eqs. (1) and (2) may be 
derived by a variational procedure if the validity of Eqs. (3) and 
(4) is assumed. Since Biot applied the variation to the displace- 
ments “4; and s;, his procedure appears to correspond to Green’s 
method. 

Thus, the question may be asked whether it might not be 
feasible to establish, for thermoelasticity, an associated principle 
of complementary energy which would then correspond to 
Castigliano’s method in elasticity. The availability of such an 
extremum principle could be of considerable significance in the 
development of approximate methods of thermoelastic structural 
analysis. 

In what follows, it is shown that it is, in fact, possible to formu- 
late such a principle complementary to Biot’s procedure. For the 
sake of brevity, the principle will not be presented here in its most 
general form. This applies in particular to the restriction in the 
type of boundary conditions. 

We consider the functional 


I= [W + D]dV — B0;) — (5) 


and set its first variation 6J equal to zero by varying 7,; and @. 
D is Biot’s dissipation function, which can be expressed as a 
quadratic form of the temperature gradient 9,; using Eq. (2), 


such that this equation may also be written in the form 
Do, +s; = 0 (2’) 
In Eq. (5) 6; stands for 6n;, where n; is the outward unit vector 
normal to the boundary A. F; is the boundary traction in the 
isothermal case (@ = 0). The 2nd integral in Eq. (5) extends over 
the surface A bounding a volume V. It can be transformed by 
means of the relations 


* This investigation was supported by the USAF, through OSR of the 
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— 30,dA = [ui(7ij.; — BO.) + 


J 56dA = [s:0,; + s,,0|dV 7 


It is assumed for simplicity that displacements wu; and s, are pre 
scribed on A. 
The variation 6/ takes then the form 


alls — — + (Woo + 4 
Bu;,:)60 — — B60.) + + si) (g 


The last 2 terms of the integrand vanish because of Eqs. (1) an¢ 
(2), and the first 2 terms yield, as required, Eqs. (8) and (4). 

The complementary energy principle associated with Biot’ 
formulation is thus established. It states, in words, that among 
all possible states of stress and temperature distribution which 
satisfy the equations of equilibrium and heat conduction in th; 
interior of the body and the prescribed boundary conditions oy 
its surface, the actually occurring states relate the stresses and 
the temperature to the solid and entropy displacements in such g 
manner as to make J stationary. 

It was pointed out above that Biot’s variational procedure 
appears to correspond to Green’s method. It will be shown by the 
present author elsewhere that this procedure is actually a mixed 
one. This is due to the fact that the equation of heat conduction 
turns out to be not an equilibrium-type equation but a relation 
of the force-displacement type. The force in this case is the 
quantity conjugate to the entropy displacement. 

It will be shown further that it is possible to establish, for 
thermoelasticity, a more general variational theorem on_ both 
displacements and stresses, as was done by Reissner? for isother- 
mal elasticity (6 = 0), which yields Eqs. (1) through (4) by means 
of a variational procedure. It is possible to show, moreover, as 
was done by Reissner® again for isothermal] elasticity, that in the 
variational theorem for displacements (Green’s method) one is 
concerned with a minimum problem while in the variational 
theorem for stresses (Castigliano’s method) one is concerned with 
a maximum problem. 

In conclusion it may be mentioned that extensions of the fore- 
going theorems to anisotropic bodies and to systems subjected 
to body forces are possible. Moreover, all of the above applies 
not only to the thermoelasticity but also to elasticity of porous 
materials because of the analogy established by Biot.! | Further- 
more, the new variational procedures may prove useful also in 
the particular case of purely thermal problems, if Biot’s* approach 
to heat flow analysis is adopted. 


REFERENCES 

! Biot, M. A., Thermoelasticity and Irreversible Thermodynamics, J. Appl 
Phys., Vol. 27, No. 3, pp. 240-253, March, 1956 

2 Reissner, E., On a Variational Theorem in Elasticity, J. Math. and Phys 
Vol. 29, pp. 90-95, 1950. 

3 Reissner, E., On Variational Principles in Elasticity, Symposium on the 
Calculus of Variations and Its Applications, Am. Math. Soc., April, 1956 

4 Biot, M. A., New Methods in Heat Flow Analysts With Application | 
Flight Structures, Journal of the Aeronautical Sciences, Vol. 24, No. 12, pp 
857-873, December, 1957. 


Further Numerical Data on Blunt Bodies 


H. Melvin Lieberstein 

Physics Section, Convair, A Division of General Dynamics 
Corporation, San Diego, Calif.* 

June 2, 1958 


A‘ ADVANTAGE of computing hypersonic flows over wind- 
tunnel measurement is the ease with which very high Mach 


* The numerical work on this note was completed while the author was 
at The Ramo-Wooldridge Corporation. 


Number 
Number 
ence 1 in 
define 
from th 
tioned Ww 
put no d 
we wish | 
Furth 
ticular, 
theory 0 
for four « 


P mea 


Case 1 


spherical 
erence ; 
Case 2 
1.4: bod 
(174. 


Case 3 


= 
M 
| 4,7 = 
- 
01074! 
1662: 
ae 022165 
P2257" 
4522: 
£5196] 
+ 1415] 
02355: 
ae 
03431 
3821: 
= 
x 
-7089 


are pre 


(1) and 
| (4). 
Biot's 
among 
which 
in th 
Ons on 
es and 


such g 


cedure 
by the 
mixed 
uction 
lation 
is the 


h, for 
both 
ther. 
neans 
eT, as 
in the 
me is 
ional 
with 


fore- 
ected 
»plies 
yrous 
ther- 
in 
oach 


Appl 
hys 
ithe 


PP 


READERS’ 


Number cases can be treated. Of course, approach-flow Mach 
Number is an input parameter and the routine described in refer- 
ence 1 includes the possibility of computing flows with an effective 
,, defined by h = |y/(y — 1)] (p/p), in the shock layer different 
from the approach flow—i.e., for y # yi. The routine men- 
tioned was tailored for close approximation of dissociation effects, 
put no data were quoted for dissociated flows in reference 1, and 
we wish to supply such data here. 

Further, one is interested in current design criteria, and, in par- 
ticular, we wish to test the merit of the modified Newtonian 
theory of Prof. L. Lees. We therefore submit the following data 


for four Cases. 


r, 2) = body coordinates 
4= arctan v/a = inclination of body to approach-flow normal 
M = Mach Number at point (r, z) on body 
P-P; ratio of pressure coefficient at a given point on body 
to pressure coefficient at stagnation 


Case 1—Shock (s + 20)? = (M? — 1) (AM? + r?), M = 5.8, 
y = 7 = 1.4; body (see Fig. 4 of reference 1) is approximately 
spherical shape for which data are given in detail in the ref- 
erence; 4k = .1585. 

Case 2—Shock + 20)? = 0.8 (0.8 + r?), = 5.8, = 
14; body (see Fig. 5 of reference 1) is a rounded cone; «xh = 
0.174. 

Case 3—Shock (s + 29)? = (M? — 1) + 7r?), M = 20, = 
14, y = 1.17 (intended to approximate the gas law from real-gas 


TABLE 1 


Case 1 
10747 11972 18927 21517 
13128 +293 0235 
o2216' -14690 39036 03124 28979 
«24048 15341 42317 3370 2813 
228577 708C 50192 3985 
, 019241 »5835 4601 
2351 220198 61579 4836 7689 
«42779 24671 74671 6835 
045222 226312 278791 26067 26559 
251961 31386 9C761 25767 
568 2356 i. 2738 25119 
TABLE 2 
Case 2 
+04278 11026 -0799 20483 09954 
1/151 012479 2575 02255 9538 
023554, 015325 4086 3569 
30722 «18362 25072 4402 +8352 
34316 20139 25508 4763 
38213 022242 5938 05114 27229 
0/3428 225320 26/47 25519 27504 
0515 26032 «7076 
06123 23769 27699 26493 26681 
©7089 24530 28208 26779 26340 
1d 84 785 2512 


FORUM 


TABLE 3 


Case 3 


x. yA M ‘ (radians 
-02819 005455 | 20285 
07838 05768 | 0922 «9900 
13238 «2224 | -9716 
1684, .07090 2822 ,1979 +9547 
|.21546 08175 359% | 29275 
754 08784 23949 2760 09132 
227093 C9811 23113 
.1865C 74,29 | «497 .728 
237768 | 1.105 
| 

TABLE 4 
Case 4 
,067c3 £2049 .0620 
i 


tables for density at 100,000 ft. altitude ); body is nearly spherical; 
kh = 0.0633. 

Case 4—Shock (2 + 20)? = (MW? — 1)(.M? + 7r?), VW = 1.2, 7 = 
y1 = 1.4; only two points on the body are given. 

It will be seen that the modified Newtonian law, C = cos? @ ts 
for Case 1, about 3 per cent high; Case 2, about 3 per cent low; 
Case 3, 6-7 per cent high; and Case 4, over 20 per cent high in the 


subsonic region. 
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O* PAGE 319, sentence after Eq. (2.6), ‘‘In terms of the polar 
coordinates, Eq. (2.5) becomes,” instead of Eq. (2.4). 

On page 321, Eq. (3.10), the differential operator should be 
[p(Op/Op) | instead of [p(Op/p) J. 

On page 322, Eq. (3.19), the factors appearing in the numer- 
ator and denominator should be 

(1 — ae~'*) ge"? and (1 + ae~*%) 

instead of 


(1 — ae'®) te’? and (1 + get? 


respectively. 


* Journal of the Aeronautical Sciences, Vol. 25, No. 5, pp. 317-323. May 
1958. 
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Similar Solution of a Laminar Boundary Layer Integrating, 
With Rotational Free Stream n+12B —raA2—a-¢ 
=e fr! e dy 
Franklin D. Hains* = (1! 
Aerodynamicist, Research Division, Bell Aircraft Corporation, fi = (a + 12B)e™ fe de dy 
Buffalo, 32 
June 6, 1958 where 
q= 12B(4B? — E), r = 1/72 ! 
— BOUNDARY-LAYER EQUATION for flow over a curved or t= E — 6B, & =(n?/24) + Bn + E\ (12 26 
plane surface is By repeated partial integrations, an asymptotic expansion js o} 
— = Pe + rbyyy (1) tained for so that the first two terms of Eq. (6) become 24 
where y is the stream function, p the density, p the pressure, and f=fitfe =&t + f&-* — [(n/4) + 3BlE-! + 
vy the kinematic viscosity. x and y are the coordinates measured [39(y + 12B)? — 2E + 12B2|(&-*/6)—...} X 20 
along the surface and normal to the surface, respectively. exp [(—7°/72) — Bn? — En] (2 
A solution will be given which satisfies the boundary conditions een * 
The power series development near 7 = 0 is 16 
y 0: u 7 0, Uy a (2) f (An?/2!) (Cn3/3!) (A2y4/2!4!) 
in such a way that for x > 0 the slope of the streamlines is every- [A2 — (4/3)AC](m5/5!) + [(7/8)A* — C2? + 8AC}(n5/6!) ~ 
where finite. In order to avoid a controversy similar to that be- [(183/30)A2C — (21/5)A* + 3C?2](n7/7!) + . (14 . 12 
ree arti « q2 > 2 > , iti 
tween Glauert! and Li? over the correct boundary condition at where C = (E/12) — (B?/2) = 


7» = ©, no attempt is made to describe the ‘‘free stream’’ in the 
region x < 0 ahead of the body. A typical streamline and 
velocity profile for the accelerating flow is shown in Fig. 1. 


The power series Eq. (14) and the asymptotic expansion Eq. (1) 
are joined at a suitable point 7 = 7* where both series are valid 
The three constants A, B, and E are found by making the value 4 


Introducing 
phe, = * 4 » hj » 
= = (3) of f , and f”’ of both series equal aty = The higher deriys. 
si tives are then automatically satisfied. 
Eq. 1 is transformed into 0 
REFERENCES 
(77/2) — = (9/2) (4) 
' Glauert, M. B., The Boundary Layer in Simple Shear Flow Past a Fi, 
and Eq. (2) becomes Plate, Journal of the Aeronautical Sciences, Readers’ Forum, Vol. 24, No. | 
9 pp. 848, 849, November, 1957. 
7 0: f f 0, at f = 1/12 (5) 2 Li, Ting Yi, Author’s Reply, Journal of the Aeronautical Sciences, Read 
The solution to Eq. (4) is obtained in the form of a power series ers’ Forum, Vol. 24, No. 11, p. 649, 850, November, 1957. Fic. 1. 
near » = 0 and an asymptotic expansion for large values of 7. 
The asymptotic expansion is taken to be * 
ow. fer Reh > (6) 
The first term, suggested by Eq. (5) is ‘ 
° Some Experiments on Effects of Yaw on 
= 2/¢ . 
Si = (n?/24) + Bn + E (7) Boundary-Layer Development in Supersonic 16( 
* 
Because f( ©) = fi( ©), Eq. (4) becomes Cone Flow 
1/2)[B? — (E£/6) — f’?] =0 8) : 
f ff + A? = fl ( D. G. DeCoursin and W. S. Bradfield 14( 
As a first approximation, Rosemount Aeronautical pomerie, University of —_ 
. Rosemount, Minn.; and Senior Staff Scientist, Scientific Researct 
Pet Fen. Fen th’. F" eh (9) Laboratory, Convair, A Division of General Dynamics Corp 12¢ 
and Eq. (8) becomes San Diego, Calif.; Respectively 
June 16, 1958 
fal” + [(n2/24) + Bn + Elfe” = 0 (10) 
- . ory NOTE is a summary of a study of the effect of small angles 
* Presently, Aerodynamicist, Boeing Scientific Research Laboratory, of yaw on boundary-layer development on a cone! including 
turbulent boundary-layer data previously obtained in connection 


4 with a more general experimental program.” * Details of the Yn 
experiments are given in references 1 and 2. = 6 

Ay In Fig. 1 the effect of yaw on the laminar velocity boundar) = 

layer on the windward side of a cone is presented. It can be ed 
seen that at the meridianal plane as the angle of yaw increases 4 


the boundary layer thins and the velocity gradient becomes 


steeper. 

On the leeward meridianal plane (Fig. 2), however, the effect is 2! 
markedly to increase the thickness of the boundary layer and 
even to change the apparent character of the boundary layet | 


from laminar to tutbulent. This effect has occurred already 4 
the angle of attack is changed from zero to 1.2° as shown in the | 


figure. 
FREESTREAM By plotting the shape parameter for these apparently turbulen! 
boundary-layer profiles for comparison with other profiles whic! Fic. 2 


were taken from the normally turbulent boundary layer at zer 
vaw, it was found’ that the shape parameters thus determined 
fell together with normal turbulent boundary-layer shape param 


iS BOUNDARY LAYER eters. Therefore, one inclines to the belief that normal turbulent | 


* The experiments presented here were carried out under Contract Num 
Fic. 1. Boundary layer with rotational free stream. ber AF18(600)-384, sponsored by the USAF, Office of Scientific Research 
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Fic. 1. Effect of yaw on laminar velocity profiles on the wind- U, 
ward surface of a 15° cone at Ma = 2.73. Fic. 3(a). Effect of variation of Reynolds Number at constant 
yaw. 
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Num Fic. 3(b). Effect of variation in yaw at constant Reynolds 
Number. 
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Fic. 4. Effect of yaw on local skin friction coefficient for laminar 
flow over a 15° cone at Wa = 2.73. 


mixing characterizes the flow along the upper meridianal plane, 
at least at these small angles of yaw. It is of interest to note 
that the Reynolds Number R, = 0.27 X 10° for these experi- 
ments is far below that at which natural transition to turbulent 
flow occurs for smooth unyawed cones at this Mach Number. 
Obviously, it is the crossflow which triggers transition. Since 
the crossflow is an important factor here, one might expect Rg to 
be a more meaningful correlation parameter in conjunction with 
some crossflow pressure gradient parameters than Reynolds 
Number based on x. 

The results of Fig. 2 are replotted on Fig. 3(b) in the y/@ di- 
mensionless representation. It is clear that the profiles are simi- 
lar when so plotted. 
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Shown for comparison [Fig. 3(a)] are profiles taken at constagf 
yaw but for a wide range of Reynolds Number.? The Mag 
Number in this case was 1, = 3.2. Again, the profiles are simj 
lar despite the fact that measurements on both the windwagg 
and leeward sides are included. 

Both sets of data are fitted by power profiles and are compared 
with the 1/7 profile chosen as characteristic of unyawed data jg 
this range. For the data of Fig. 3(b), a 1/10 power profile seems 
suitable whereas for the data of Fig. 3(a), a 1/8 profile works vegy 
well. In the case of the 1/10 power profile, the crossflow pressu 
gradient is greater than for the data of Fig. 3(a) because the 
measuring station was closer to the cone tip. 

Fig. 4 shows the experimental variation of laminar skin frictiog 
coefficients on the windward side with vaw at 1/,, = 2.73 on the 
15° total angle cone. It can be seen from these results that the 
skin friction coefficient varies rapidly for 0 < a@ < 1 and thaf 
further increases in yaw have a less strong effect on skin frictiog 
coefficients. A comparison of these experimental results with 
theoretical analyses is shown. 

When the angle of yaw of the cone increased to approximate 
one half the nose angle, separation of the boundary layer from 
the upper surface of the body starting from the tip occurred.! 
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Boundary-Layer Equation for Axially Symmetric Flow . . . 


(Continued from page 634) 


is due to different methods of calculation. If our 
method of calculation is applied to Tomotika’s data, 
the corresponding curve of friction obtained coincides 
with ours. 

Up to 6 = 50°, we used the values of a corresponding 
to Hartree’s solution; therefore, the curve of friction 
could not have lower values. Fage‘ calibrated his ap- 
pliance to measure friction by using Tomotika’s calcu- 
lations. This should account, at least partially, for the 
difference between Fage’s results and ours. 
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